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INTRODUCTION 

While  manual  palpation  remains  the  first  diagnostic  line  of  defense  against  breast 
cancer[l-4],  it  is  unfortunately  limited  to  relatively  large  and  superficial  lesions.  The 
central  hypothesis  of  this  work  is  that  remote  measurement  of  elasticity  in  breast  tissues 
is  possible  and  provides  unique  information  which  could  increase  detection  and/or 
characterization  of  potentially  malignant  masses  not  manually  accessible.  Our 
preliminary  studies  suggest  that  the  proposed  methods  are  capable  of  precisely  measuring 
internal  deformation  and  strain  in  three  dimensions.  These  data  are  required  to 
reconstruct  the  elasticity  distribution  within  the  object.  Consequently,  technologies 
developed  within  the  scope  of  this  project  may  have  significant  diagnostic  value  in 
detection  and  management  of  breast  cancer. 

Most  elastography  to  date  has  utilized  ultrasound  [5-14],  although  MRI  has  also 
been  used  [15-22].  Usually  an  external  static  or  dynamic  deformation  is  applied  while 
the  resultant  displacement  or  propagating  shear  wave  is  documented  by  imaging  devices. 
In  general,  to  reconstruct  the  tissue-specific  property  of  Young’s  modulus  in  complex 
systems  such  as  the  breast,  the  3-D  displacement  vector  must  be  measured  over  a  3-D 
volume.  In  this  project,  we  are  developing  techniques  that  measure  the  3-D  displacement 
vector  over  any  volume  in  the  object  at  high  spatial  resolution.  These  data  will  be 
processed  to  produce  3-D  strain  images  for  submission  to  a  3D  elasticity  reconstruction 
routine  to  map  the  relative  Young’s  modulus  within  the  volume. 

ACCOMPLISHMENTS  IN  RELATION  TO  STATEMENT  OF  WORK 

The  goal  of  this  research  program  is  to  develop  a  sensitive  diagnostic  technique 
based  on  quantitative  elasticity  imaging  permitting  surrogate  palpation  of  deep  lying 
breast  lesions.  To  this  end,  specific  technical  objectives/tasks  were  adopted  and 
summarized  in  the  proposal  Statement  of  Work.  Each  task  is  restated  below  along  with 
descriptions  of  relevant  Methods  and  Progress/Results. 

Technical  Objective  I:  Data  Acquisition  and  Reconstruction 
Task  Is  Month  1-3 

Construction  of  computer  controlled  hydraulic  compression  device  capable  of  producing 
an  incremental  surface  deformation  of  the  mechanical  body.  This  device  will  be 
triggered  by  NMR  imaging  system  and  have  a  simple  timing/displacement  control 
appropriate  for  phantom  studies. 

The  proposed  method  requires  that  the  imaged  object  experience  a  relatively 
minor  externally-applied  deformation  force.  Based  on  signal  optimization  simulations, 
differential  deformations  below  10%  (e.g.,  <lcm  deformation  across  a  10cm  object) 
should  be  adequate  for  generation  of  elasticity  maps.  Design  details  of  the  mechanical 
deformation  device  were  provided  in  Year  1  Progress  Report,  although  are  briefly 
summarized  here  with  the  device  schematic  is  shown  in  Figure  1.  Four  neoprene  bellows 
outside  of  the  NMR  coil  (rf-coil)  provide  both  upward  and  downward  force  to  an  acrylic 
plate  on  top  of  the  phantom.  The  top  bellows  are  on  a  common  pressure  circuit,  likewise 
the  bottom  pair  are  pressurized  together  to  yield  uniform  vertical  displacement  along  the 
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top  of  phantom.  Solid  acrylic  yolks  at  both  ends  of  the  device  insure  that  only  the  top 
plate  moves  in  response  to  pressurization.  Pneumatic  elements  are  driven  by  air-pressure 
solenoid  values  that  are,  in  turn,  controlled  by  a  transistor-transitor-logic  circuit  triggered 
by  the  NMR  imaging  sequence. _ 

Pressurize  for 
MAX  deformation 


Pressurize  for 
MIN  deformation 


Figure  1.  Side-view  schematic  of  pneumatic  deformation  system. 


Phase-encode  MRI  methods,  as  employed  here,  require  that  data  are  acquired  in 
many  segments  over  an  extended  period  of  time.  Consequently  the  mechanical 
deformation  must  be  highly  reproducible  for  each  segment.  Mechanical  stability  tests  of 
the  device  were  performed  using  a  slightly  modified  version  of  the  2D-displacement 
encoding  MRI  sequence.  The  upshot  of  these  tests  was  that  the  system  is  extremely 
stable  and  reproducible  in  application  of  a  5mm  differential  deformation  (mixing  time, 
TM=350msec).  The  observed  phase  instabilities  were  barely  above  that  allowed  by  NMR 
signal  limits.  These  results  suggest  mechanical  reproducibility  within  <20microns.  One 
significant  finding  of  these  stability  tests  was  that  the  mechanical  response  of  the 
deformation  device  was  significantly  faster  (~  150ms)  than  the  “settling  time”  of  the 
tissue-mimicking  phantoms  (~350ms).  That  is,  internal  mechanical  reflected  waves 
persist  beyond  point  the  acrylic  drive  plates  stop.  We  anticipate  that  breast  tissue  will 
have  comparably  long  settling  times,  thus  we  will  continue  to  use  long  mixing  times 
(~350ms)  in  our  experiments  despite  the  SNR  gains  of  shorter  times.  Immunity  to  long 
persistence  mechanical  waves  was  the  main  motivation  to  use  the  stimulated-echo 
approach  as  is  done  here. 


Task  2:  Month  1-6 
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Development  of  various  phantoms  to  test  NMR  displacement  and  strain  data  acquisition 
and  3-D  elasticity  reconstruction  algorithms. 

Fabrication  of  phantoms  to  model  the  mechanical  and  NMR  properties  of  human 
tissue  is  an  important  step  in  development  and  verification  of  optimized  NMR 
displacement-encoded  volume  simulated-echo  pulse  sequence.  Moreover,  simple 
geometries  of  phantoms  allow  characterization  of  the  elasticity  reconstruction  algorithms. 
The  phantom  materials  should  closely  resemble  relevant  tissue  properties.  In  addition, 
these  materials  should  be  stable  (i.e.  months)  and  permit  the  non-destructive  embedding 
of  lesion-equivalent  targets  within  the  phantom. 

We  have  previously  developed  tissue-mimicking  phantoms  to  test  NMR  elasticity 
imaging  device.  These  phantoms  were  made  of  two  materials:  a)  polymer  produced  by 
M-F  Manufacturing  Co.,  Inc.  (Fort  Worth,  TX),  and  b)  gelatin  (Sigma-Aldrich  Co.,  3050 
Spruce  Street,  St.  Louis,  MO  63103  USA).  The  first  material,  plastisol,  consists  of  a 
liquid  plastic  combined  with  either  softener  (plasticizer)  or  hardener.  By  varying  the 
proportion  of  these  two  components,  it  is  possible  to  produce  composite  models  of 
desired  elasticity.  The  raw  composite  materials  were  stirred  and  heated  to  approximately 
170°C.  At  that  temperature  the  mixture  was  poured  into  molds  producing  a  tissue- 
mimicking  time-stable  phantom  of  desired  shape  and  elasticity  distribution.  However, 
two  major  deficiencies  remain.  First,  plastisol  does  not  have  desired  tissue  equivalent 
NMR  properties.  Second,  due  to  high  temperature  rise  during  phantom  preparation,  it  is 
impossible  to  produce  tissue-containing  phantoms  using  plastisol. 

Several  different  materials  were  further  considered  for  NMR  elasticity  phantoms. 
These  materials  include  Semicosil  921,  Semicosil  905  and  Silgel  612  (Wacker  Silicones 
Co.,  Adrian,  MI  49221),  and  Rhodorsil  RTV  163  (Rhone-Poulenc,  France).  From  all 
tested  materials,  the  Semocosil  921  silicone  gel  appeared  to  mimic  the  tissue  properties 
the  best.  The  elasticity  of  silicone  gel  can  be  simply  controlled  by  mixing  ratio  of  two 
components.  The  composite  material  is  then  poured  into  the  mold  and  cured  at  the  room 
temperature  enabling  fabrication  of  tissue  containing  phantoms.  These  silicone  gels  have 
shown  high  SNR  spin-echo  signal  over  the  range  of  hardness  variations.  In  addition,  the 
mechanical  properties  of  silicone  did  not  change  over  60  days,  as  was  measured  by 
Instron-type  mechanical  system.  Therefore,  it  is  possible  to  fabricate  tissue-mimicking 
phantoms  of  desired  shape  and  elasticity  distribution  using  Semicosil  921  silicone  gel. 

Based  on  favorable  NMR  and  mechanical  properties,  phantoms  were  fabricated 
using  the  Semicosil  materials.  These  phantoms  were  designed  with  simple  geometries  to 
assess  spatial  resolution  and  accuracy  of  3D  strain  imaging  3D  elasticity  reconstruction 
algorithms,  and  include  crossed  bars  of  hard  material  (6x  Young’s  modulus)  with  variable 
size  and  separation  (for  resolution)  and  simple  spherical  inclusions  for  comparison  with 
theoretical  results. 

Task  3:  Month  1-12 

Development  of  the  NMR  displacement-encoded,  volume  stimulated-echo  pulse  sequence 
optimized  for  displacement/strain  sensitivity  and  SNR. 

Two,  3D  displacement-encode,  volume  stimulated-echo  pulse  sequences  have 
been  developed.  These  are  shown  in  Figure  2(a)  using  "classic"  3D  spatial  encoding,  and 
(b)  a  fast-spin-echo  (FSE)  extension  to  more  efficiently  encode  the  3rd  dimension.  The 
classic  3D  sequence,  is  simplest  and  implement  but  unfortunately  it  is  clinically 
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impractical  due  to  very  long  scantimes.  For  example,  acquisition  of  volumetric  data  at 
256x128  matrix  resolution  over  16  slices  with  full  3-axis  displacement  sensitivity 
typically  requires  l-2hours  (if  TR=lsec;  TRx  128x1 6x3= 1.7hrs).  Incorporation  of  the 
FSE  echo  train  permits  a  scantime  reduction  by  a  factor  of  1/4,  1/8,  to  1/16  dependent  on 
echo  train  length.  To  date,  we  have  implemented  an  8-echo  FSE  sequence  where  the  FSE 
phase-encoding  is  applied  along  the  3rd  dimension  (i.e.  slice  direction).  Ideally,  all  slice 
encode  steps  are  acquired  in  one-shot  which  greatly  reduces  scan  time  (to  13min). 
Unfortunately,  it  appears  that  poor  gradient  hardware  performance  is  introducing  eddy 
current  effects  that  confound  3rd  axis  encoding  on  our  2T  MRI  system.  Phase  shifts  from 
this  artifact  derail  strain  measurement  based  on  phase.  We  are  in  the  process  of 
evaluating  means  to  apply  phase  corrections  to  reduce  imperfections  in  the  FSE 
component  of  the  sequence. 


Figure  2.  Volumetric  displacement-encoding  stimulated  echo  sequences  using  (a)  classic 
3D  encoding  and  (b)  fast-spin-echo  encoding  of  the  third  dimension. 
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Images  of  tissue  mimicking  phantom  containing  ramped  bars  of  harder  material  acquired 
via  the  3D  FSE  stimulated  echo  sequence  are  shown  in  Figure  3.  Spurious  intensity 
modulations  in  the  magnitude  image  (i.e.  conventional  MRI  contrast)  in  Figure  3(a) 
indicate  an  imperfect  reconstruction  along  the  3rd  dimension.  That  is,  there  is  an 
interference  of  signal  from  adjacent  slices.  Despite  this,  derivative  displacement  and 
strain  images  display  features  of  the  object  (Figure  3(b-f)),  however,  residual  phase  error 
prevent  acceptable  elasticity  reconstructions.  Deformation  and  object  symmetry  meant 
the  displacement  and  strain  in  the  3rd  dimension  (ie.  perpendicular  to  the  image)  was 
small.  This  material  was  presented  at  the  1999  International  Society  of  Magnetic 
Resonance  Medicine  (ISMRM)  and  is  included  in  the  appendix. 
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Figure  3.  Volumetric  displacement-encoding  stimulated  echo  sequence  applied  to  Semicosil  rubber 
phantom  containing  crossed  ramped  bars  of  hard  material.  Shown  are  (a)  magnitude  image, 
displacement  images  in  (b)  X  right-to-left,  (c)  Y  top-to-bottom,  and  (d)  Z  through-plane.  Strain 
images  (i.e.  spatial  derivative  of  displacement)  for  X  and  Y  displacements  are  shown  in  (e)  and  (f) 
respectively. 

Computer  simulations  were  performed  to  address  the  tradeoff  between  displacement 
sensitivity  and  NMR  signal-to-noise  (SNR)  since  optimization  of  one  parameter  comes  at 
the  expense  of  another.  Key  parameters  that  affect  SNR  and  displacement  sensitivity 
include:  displacement  gradient  amplitude  and  duration;  applied  deformation  amplitude; 
image  resolution;  mixing  time,  TM;  and  tissue  water  diffusion  properties.  Results  of  the 
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optimization  studies  suggest  greater  SNR  can  be  achieved  through  greater  surface  strain 
(i.e.  deformation  of  the  object  by  -10%)  with  a  corresponding  reduction  in  displacement 
encoding  (reduced  from  4  to  1  gauss/cm).  Increasing  the  applied  deformation  will  have 
the  added  benefit  to  make  the  applied  deformation  more  significant  relative  to 
"background"  physiologic  motions  such  as  cardiac  and  respiratory.  Yet  higher 
deformation,  however,  introduce  the  need  to  consider  non-linear  elastic  effects  (see  Task 
5).  These  optimization  studies  were  also  presented  at  the  1999  International  Society  of 
Magnetic  Resonance  Medicine.  The  abstract  is  included  in  the  Appendix. 


Task  4:  Month  1-12 

Expansion  of  previously  developed  linear  elasticity  reconstruction  algorithms  for 
volumetric  displacement  and  strain  NMR  measurements. 

Expansion  of  linear  elasticity  reconstruction  algorithms  were  focused  on  the 
following  two  aspects:  a)  development  of  boundary  detection  methods  to  identify  the 
regions  of  uniform  elasticity  (Young’s  modulus)  distribution,  and  b)  development  of 
volumetric  linear  elasticity  reconstruction  algorithms  for  optimized  NMR  displacement- 
encoded,  volume  simulated-echo  pulse  sequence. 

The  boundary  detection  algorithm  was  developed  for  specifically  NMR  elasticity 
imaging.  This  algorithm  is  based  on  the  stress  continuity  condition  applicable  for  soft 
tissue  deformations  used  in  these  studies.  The  relevant  boundary  definition  methods  were 
detailed  in  the  Year  1  Progress  Report. 

The  following  relates  to  extensions  to  a  3D  elasticity  reconstruction  algorithm 
wherein  complete  three-dimensional  strain  data  are  required  to  solve  for  a  general,  three- 
dimensional  object.  Previously,  only  2D  reconstruction  algorithms  have  been  applied, 
yet  these  are  anticipated  to  lead  to  errors  at  planes  where  through-plane  strain  is 
significant.  As  an  illustration.  Figure  4  shows  simulated  elasticity  distributions  at  a  slice 
just  within  (z=0.95  radius)  and  just  outside  (z=1.05  radius)  of  a  hard  spherical  inclusion. 
These  data  were  used  to  test  the  3D  elasticity  reconstruction  algorithm  and  provide 
supportive  evidence  for  the  need  of  3D  reconstruction  algorithms  in  general. 


(a)  <c)  (d)  Centra]  profiles 
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Figure  4.  Simulated  elasticity  maps  at  planes  just  outside  of  (top  row)  and  just  within  of  (bottom 
row)  a  spherical  hard  inclusion.  The  3D  reconstruction,  k3,  (images  on  right)  was  a  more  faithful 
depiction  of  truth,  k,  (images  on  left),  relative  to  the  2D  reconstruction,  k2,  (images  in  middle). 
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Task  5:  Month  10-26 

Development  of  the  nonlinear  elasticity  reconstruction  model  capable  to  process  high- 
strain  NMR  images. 

The  elasticity  reconstruction  problem  can  be  posed  in  number  of  ways  depending 
on  the  experimental  conditions  such  as  amount  of  surface  deformations,  type  of  phantom 
materials,  etc.  Previously,  we  have  developed  a  reconstruction  algorithm  based  on  linear 
theory  of  elasticity,  i.e.,  applicable  for  small  deformations  and  linear  elastic  materials. 
However,  both  mechanical  equilibrium  equations  and  the  displacement-strain  relations 
must  include  high  order  spatial  derivatives  of  the  displacement  if  the  average  deformation 
is  too  large  to  be  described  with  a  linear  model.  An  algorithm  was  developed  to 
reconstruct  the  elastic  modulus  of  soft  tissue  for  comparatively  large  surface 
deformations  required  for  optimal  SNR  and  sensitivity  strain/displacement  images. 
Numerical  methods  were  developed  to  reduce  error  propagation  in  reconstruction 
algorithms. 

As  said,  linear  elastic  theory  is  reasonable  in  the  realm  of  a  few  percent  surface 
deformation  as  was  used  in  our  initial  experiments.  Based  on  optimization  studies 
performed  for  Task  3,  however,  greater  sensitivity  and  SNR  is  available  at  higher 
deformation,  such  as  >10%.  Consequently,  in  current  experiments  we  employ  surface 
deformations  -10%.  As  predicted,  at  these  settings  there  is  sufficient  phase  contrast  of 
internal  deformations  with  reduced  signal  losses  to  diffusion  since  displacment-encoding 
gradients  are  proportionally  reduced.  Unfortunately,  image  artifact  due  to  imperfect 
hardware  performance  on  our  2T  unit  using  the  3D  fast-spin-echo  (FSE)  exceed  non¬ 
linear  effects.  Moreover,  since  the  artifact  is  unrelated  to  linearity  of  elastic  deformation, 
application  of  non-linear  algorithm  elements  is  not  a  remedy  to  their  removal.  We  have 
decided  to  forego  application  of  the  non-linear  modification  of  the  reconstruction 
algorithm  until  we  resolve  the  artifacts  in  the  3D  FSE  technique. 

Technical  Objective  II:  Phantom  Studies  on  2T  18cm  Bore  MRI  System 

Task  6:  Month  7-22 

Development  of  gel-  and  rubber-based  phantoms  with  tissue  mimicking  elastic  properties. 
Time  stable  phantoms  with  non-palpable  inclusions  of  various  shapes  and  elasticity 
contrast  positioned  at  different  locations  within  the  phantom  will  be  produced. 

During  the  first  year  of  the  project,  we  have  identified  and  tested  the  materials  for 
tissue  mimicking  phantoms.  Our  ultimate  goal  is  to  simulate  the  anatomical  and 
geometrical  features  of  the  normal  and  pathologically  transformed  breast  using  these 
materials.  The  models  of  breast  containing  single  or  multiple  inclusions  were  fabricated 
using  plastisol  material  (M-F  Manufacturing  Co.,  Fort  Worth,  Texas,  USA). 

Our  recent  studies,  however,  suggested  that  silicone  gels  could  better  simulate  the 
mechanical  and  NMR  relaxation  properties  of  the  tissue.  Initially,  silicone-gel  based 
homogeneous  phantoms  were  produced  for  mechanical  and  NMR  testing  (Semicosil  921, 
Wacker  Silicones  Co.,  Adrian,  Michigan,  USA).  These  tests  showed  the  material  had 
suitable  NMR  properties  -  like  tissue.  To  control  mechanical  properties,  the  Semicosil 
921  contains  two  components,  A  and  B,  wherein  different  ratios  of  these  components  are 
used  to  vary  the  mechanical  properties  of  the  gel.  A  tissue-mimicking  phantom  was 
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constructed  in  several  steps.  First,  background  material  was  prepared  by  thoroughly 
mixing  components  A  and  B  in  a  1:1  ratio,  and  then  pouring  the  mixture  into  a  154-mm 
by  80-mm  rectangular  mold.  The  mixture  was  degassed  and  cured  for  24  hours  at  room 
temperature  to  produce  a  22-mm  thick  layer.  Then  a  25-mm  diameter  hard  sphere  was 
prepared  from  a  1:2.5  mixture  of  A  and  B  and  was  placed  on  top  of  the  layer  in  the 
middle  of  the  mold.  Finally,  another  batch  of  background  material  (1:1  ratio)  was  poured 
into  the  mold  resulting  in  a  64-mm  by  80-mm  by  154-mm  phantom  with  a  single,  hard, 
spherical  inclusion  roughly  in  the  center.  At  the  same  time,  three  samples  of  each  batch 
were  taken  to  independently  assess  the  elasticity  contrast  between  the  inclusion  and 
surrounding  materials.  These  measurements  were  performed  using  the  force-deformation 
system  described  in  (Erkamp  \etal  1998),  and  showed  that  the  inclusion  was  four  times 
harder  than  the  background,  and  that  both  background  materials  were  elastically 
equivalent.  These  are  highly  stable  and  are  still  in  use  several  months  post-fabrication. 

Task  7:  Month  12-26 

Investigation  of  the  capabilities  of  3-D  NMR  elasticity  imaging,  i.e.,  determine 
sensitivity,  accuracy  and  resolution  of  NMR  displacement  and  strain  images,  and 
reconstructed  elasticity  (Young’s  modulus)  distribution. 

As  shown  by  simulation  presented  in  Figure  4,  elasticity  maps  reconstructed  from  3- 
D  strain  data  should  be  more  accurate  than  a  2-D  elasticity  reconstructions.  This 
should  be  particularly  true  near  the  edges  of  an  object  where  the  plane-strain-state  is 
violated;  that  is  where  the  strain  along  the  “through-plane”  direction  in  non- 
negligible.  This  concept  was  tested  on  silicone-gel  based  phantoms  developed  in 
Task  6  using  the  volume  stimulated-echo  pulse  sequence  (Task  3).  A  comparison  of 
the  actual  “2D”  reconstruction  and  “3D”  reconstructed  elasticity  maps  for  2  of  32 


(a)S1(x1,x2)  (t>)  (c)  (d)  Central  profiles 


10  20  30 


Figure  5.  Reconstructed  elasticity  maps  at  planes  just  outside  of  (top  row)  and  just  within 
(bottom  row)  a  spherical  hard  inclusion  in  a  silicone-gel  phantom.  These  results  are  consistent 
with  the  simulations  (Figure  4)  that  show  the  3D  reconstruction,  k3,  (images  on  right)  more 
faithfully  depicts  the  size  of  the  inclusion. 
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planes  through  the  phantom  data  are  shown  in  Figure  5.  The  planes  selected  through 
the  phantom  and  general  format  of  Figure  5  were  chosen  for  direct  comparison 
between  the  simulation  (Figure  4)  and  experiment  (Figure  5).  Ideally,  there  would 
be  no  inclusion  visible  in  the  plane  outside  of  the  inclusion  (ie.  reconstructions  in  the 
top  row  of  Figure  5  should  be  uniform).  Violation  of  the  non-negligible  strain 
remote  from  the  object,  however,  results  in  inaccuracies  in  2D  reconstructions. 

These  errors  are  clearly  reduced  in  the  3D  reconstruction  which  confirms  simulation 
results. 


Technical  Objective  IH:  Translation  to  1.5T  Human  MRI  System 
Task  8:  Month  20-25 

Design  and  construction  of  the  compression  device  compatible  with  human  breast 
NMR  imaging  system  and  capable  to  produce  wide  range  of  surface  deformations. 
Development  of  sophisticated  circuits  to  control  surface  deformations  and  time 
synchronization  with  human  MRI  system. 

Translation  to  the  1.5T  unit  system  has  begun.  We  have  a  compressor  pneumatic- 
drive  system  in  hand  that  is  suitable  for  operation  on  our  research  1.5T  human  MRI 
unit.  Existing  phantoms  will  be  used  for  testing  on  the  1.5T  unit. 

Task  9:  Month  24-34 

Development  of  appropriate  for  clinical  studies  3-D  displacement  encoded,  volume 
stimulated-echo  pulse  sequence  on  human  MRI  system  with  data  acquisition  time 
within  or  comparable  to  regular  MRI  examination. 

Two-  and  three-dimensional  fast-spin-echo  (FSE)  are  being  modified  for 
displacement-sensitive,  stimulated-echo  acquisitions.  We  expect  to  initial  sequence 
testing  within  a  month.  When  a  32-echo-train-length  sequence  is  realized,  without 
significant  system-related  artifact,  acceptable  scan  times  should  be  achieved. 

Task  10:  Month  22-34 

Development  of  time  efficient  elasticity  reconstruction  algorithms  more  suitable  for 
clinical  applications. 

Our  current  3D  elasticity  reconstruction  is  extremely  long  (many  hours).  Since  the 
reconstruction  is  a  numerical  solution  to  invert  high-order  differential  equations, 
reconstruction  time  scales  with  the  acquired  quantity  of  slices  and  the  in-plane 
resolution.  Previously,  volumetric  datasets,  as  generated  in  this  project  did  not  exist, 
thus  not  all  practical  issues  related  to  3D  elasticity  reconstruction  were  known.  To 
date,  most  programming  efforts  have  been  directed  to  management  of  the  acquired 
six-dimensional  datasets  (i.e.  position  in  x,y,z,  and  displacement-sensitive  ux,uy,uz) 
which  expand  to  and  a  “12- vector”  per  pixel  (position,  and  9-elements  of  the  strain 
tensor).  The  12-dimensional  data  is  only  the  input  to  the  reconstruction. 
Unfortunately,  given  these  management  issues,  to  date  there  has  not  been  significant 
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progress  toward  improving  reconstruction  efficiency.  Greater  interest  is  being  paid  to 
validation  of  the  reconstructions  rather  that  reconstruction  speed. 


Task  11:  Month  20-30 

Estimate  the  influence  of  the  cardiac  and  respiratory  motion  to  elasticity  images  and 
develop  the  approaches  to  reduce  these  artifacts. 

Resolution  of  systematic  artifacts  that  originate  from  imperfect  hardware 
performance  are  more  pressing  than  estimation  of  patient-originated  artifact.  As 
such,  action  on  Task  11  are  postponed  until  we  have  an  operational  sequence  on  the 
1.5T  unit.  We  anticipate  significantly  greater  performance  on  the  1.5T  unit  which 
routinely  provides  32-echo-train-lenegth  images  without  apparent  artifact. 


Task  12:  Month  30-36 

Validation  of  clinical  NMR  data  acquisition  and  elasticity  reconstruction  methods  on 
breast  tissue  mimicking  phantoms. 


CONCLUSIONS 

Tasks  to  design,  fabricate  and  refine  the  “deformation  device”,  essential  for 
phantom  studies,  is  complete.  The  MRI-compatible  hardware,  pneumatic  components, 
and  control  circuitry  are  now  fully  operational.  This  device  provides  excellent  control 
and  stability  in  repetitive  deformations  of  an  tissue-mimicking  objects.  Significant  effort 
was  directed  to  develop  suitable  phantom  materials.  Requirements  include:  temporally 
stable  (over  months-years);  adjustable  mechanical  properties  to  match  a  range  of  tissues; 
moldable  to  simple  geometries;  tissue-like  NMR  properties.  The  Semocosil  921  silicone 
gel  generally  meets  these  requirements  (with  the  exception  of  tissue-like  diffusion)  and 
will  be  used  in  subsequent  phantom  studies.  Three-dimensional,  stimulated-echo 
acquisition  sequences  that  have  sensitivity  to  3 -dimensional  displacement  have  been 
written  and  applied  to  gather  3D  strain  data  on  simple  objects.  These  data  have  been 
input  in  newly-designed  3D  elasticity  reconstruction  routines  to  yield,  to  the  best  of  our 
knowledge,  the  first  elasticity  reconstruction  based  on  volumetric  internal  spatial/strain 
data.  Methods  to  reduce  data  acquisition  time,  via  stimulated-echo  spatial  encoding  have 
not  yet  been  satisfactory  in  terms  of  artifact  control.  We  believe  these  artifacts  are  due  to 
the  particular  MRI  unit  used,  and  have  thus  begun  to  move  experiments  to  the  1.5  T 
human  MRI  unit  which  is  know  to  have  better  performance. 
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Elasticity  Reconstructive  Imaging  by  Means  of 
Stimulated  Echo  MRI 

Thomas  L.  Chenevert,  Andrei  R.  Skovoroda,  Matthew  O’Donnell, 
Stanislav  Y.  Emelianov 


A  method  is  introduced  to  measure  internal  mechanical  dis¬ 
placement  and  strain  by  means  of  MRI.  Such  measurements 
are  needed  to  reconstruct  an  image  of  the  elastic  Young’s 
modulus.  A  stimulated  echo  acquisition  sequence  with  addi¬ 
tional  gradient  pulses  encodes  internal  displacements  in  re¬ 
sponse  to  an  externally  applied  differential  deformation.  The 
sequence  provides  an  accurate  measure  of  static  displace¬ 
ment  by  limiting  the  mechanical  transitions  to  the  mixing 
period  of  the  simulated  echo.  Elasticity  reconstruction  in¬ 
volves  definition  of  a  region  of  interest  having  uniform 
Young’s  modulus  along  its  boundary  and  subsequent  solution 
of  the  discretized  elasticity  equilibrium  equations.  Data  acqui¬ 
sition  and  reconstruction  were  performed  on  a  urethane  rub¬ 
ber  phantom  of  known  elastic  properties  and  an  ex  vivo  ca¬ 
nine  kidney  phantom  using  <2%  differential  deformation. 
Regional  elastic  properties  are  well  represented  on  Young’s 
modulus  images.  The  long-term  objective  of  this  work  is  to 
provide  a  means  for  remote  palpation  and  elasticity  quantita¬ 
tion  in  deep  tissues  otherwise  inaccessible  to  manual  palpa¬ 
tion. 

Key  words:  elastic  Young’s  modulus;  magnetic  resonance 
imaging;  elastography;  strain  imaging. 

INTRODUCTION 

It  is  well  known  that  tissue  elastic  properties  may  be 
altered  by  tumors.  Young’s  elastic  moduli  may  differ  by 
orders  of  magnitude  in  soft  tissues  in  various  physiologic 
states  (1,  2).  This  finding  is  the  physical  basis  behind 
manual  palpation  used  to  detect  “hard”  masses  (3,  4). 
Indeed,  physical  examination  is  the  first  diagnostic  line 
of  defense  against  breast  cancer,  because  nodule  hard¬ 
ness  raises  suspicion  of  malignancy.  Detection  of  a  new 
breast  mass  by  physical  examination  is  often  sufficient 
for  surgical  excisional  biopsy,  even  when  not  corrobo¬ 
rated  by  other  diagnostic  tests.  Manual  palpation  of  the 
prostate,  superficial  lymph  nodes,  and  abdominal  organs 
are  also  commonly  performed.  Unfortunately,  sensitivity 
of  palpation  is  relatively  poor  within  deep,  dense,  or 
heterogeneous  tissues.  Although  the  touch  of  a  skilled 
interpreter  is  considered  a  powerful  diagnostic  instru¬ 
ment,  most  lesions  detected  by  palpation  tend  to  be  rel¬ 
atively  large  and  superficial. 
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Scientists  are  attempting  to  electronically  extend  the 
touch  of  the  physical  examiner  using  a  variety  of  image- 
based  techniques  that  infer  tissue  elasticity.  The  essential 
element  is  measurement  of  internal  motion  and  strain  in 
tissue  structures  experiencing  mechanical  stress.  To 
date,  most  “elastography”  has  used  ultrasound  to  track 
the  relative  motion  of  targets  by  specular  reflection  (5-7), 
by  Doppler  techniques  (8-10),  by  cross-correlation  of 
raw  or  processed  acoustic  echoes  (11,  12),  or  by  tracking 
speckle  patterns  (13-15).  Usually  an  external  static  or 
dynamic  deformation  is  applied  while  internal  displace¬ 
ments  or  propagating  shear  waves  are  documented  by 
imaging. 

MRI  has  also  been  used  to  measure  internal  displace¬ 
ment  and  strain  components  of  the  heart  using  spatial 
magnetization  tagging  (16,  17)  and  phase-based  velocity 
encoding  (18).  Elasticity  reconstruction  of  an  externally 
deformed  phantom  was  demonstrated  using  magnetiza¬ 
tion  tagging,  but  this  method  has  spatial  resolution  lim¬ 
ited  by  the  tagged  grid  size  and  only  measures  2D  motion 
(19).  More  recently,  motion  phase  encoding  by  means  of 
bipolar  gradients  was  used  to  produce  two-dimensional 
(2D)  displacement  and  strain  maps  in  media  mechani¬ 
cally  driven  by  external  forces  (20-23).  Strain  and  dis¬ 
placement  maps  infer  internal  elasticity  but  are  also 
strongly  affected  by  the  applied  deformational  geometry. 
Consequently,  these  maps  do  not  uniquely  reflect  inter¬ 
nal  tissue  properties  (i.e.,  elastic  Young’s  modulus). 
Maps  of  dynamic  strain-wave  propagation,  however,  do 
allow  measurement  of  local  strain  wavelength  or  velocity 
from  which  the  local  elastic  modulus  can  be  derived  (20). 
Shear-wave  attenuation,  interference  from  standing 
waves  off  multiple  reflectors,  and  limited  resolvable 
points  over  the  shear  wavelength  are  potential  drawbacks 
of  this  approach. 

Relative  to  ultrasound,  MRI  has  an  advantage  in  overall 
resolution  and  accuracy  for  multidimensional  displace¬ 
ment  and  strain  measurement  needed  for  elasticity  re¬ 
construction.  Ultrasound  can  accurately  measure  axial 
(i.e.,  along  the  beam  axis)  motion  at  high  spatial  resolu¬ 
tion  (^millimeter),  but  lateral  displacement  is  measured 
at  much  lower  spatial  resolution  defined  by  the  depth- 
dependent  acoustic  beam  width.  The  third  dimension  is 
generally  not  even  considered  given  the  limitations  of 
ultrasound.  Consequently,  reduced  motion  dimensional¬ 
ity  and  overall  low  motion  resolution  of  the  imaging 
system  compromise  the  “elastogram.”  This  shortcoming, 
in  turn,  constrains  the  mechanical  model  used  in  elastic¬ 
ity  reconstruction.  To  date,  only  ID  motion  models  have 
been  applied  to  ultrasound-derived  elasticity  images  of 
tissues  in  vivo  (12,  24).  More  accurate  elasticity  images 
are  achieved  by  properly  controlling  external  deforma- 
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tions,  leading  to  2D  elasticity  reconstruction  within  the 
imaging  plane  (25,  26). 

In  this  work  we  present  a  method  to  spatially  encode 
internal  displacement  of  an  object  that  has  undergone  an 
externally  applied  “static”  deformation  with  subsequent 
reconstruction  into  elasticity  maps.  Unlike  dynamic 
techniques  designed  to  estimate  elasticity  from  observa¬ 
tions  related  to  strain-wave  propagation,  static  elasticity 
reconstruction  involves  estimation  of  local  strain  from 
displacement  and  numerical  solution  of  differential  elas¬ 
ticity  equilibrium  equations. 


THEORY  OF  RECONSTRUCTIVE  ELASTICITY 
IMAGING 


The  goal  of  elasticity  imaging  is  to  reconstruct  the  elastic 
modulus  of  a  desired  tissue  region  using  available  mea¬ 
surements  of  displacement  and  strain  components.  In¬ 
deed,  the  mechanical  properties  of  tissue  are  ultimately 
linked  to  the  patterns  of  internal  deformations,  but  the 
deformational  geometry  can  greatly  affect  these  patterns 
as  well.  To  uniquely  image  tissue  elasticity,  the  Young’s 
modulus  must  be  reconstructed  from  estimates  of  inter¬ 
nal  displacement  and  strain. 

In  this  paper,  the  general  approach  to  elasticity  recon¬ 
struction  was  based  on  a  model  of  linear,  elastic,  isotro¬ 
pic,  incompressible  media  (26,  27).  The  key  equations 
and  considerations  are  briefly  presented  here.  A  more 
detailed  description  of  elasticity  reconstruction  is  given 
in  an  earlier  publication  (26). 

In  linear  elasticity,  the  components  of  the  strain  (e^) 
and  stress  (cqj)  tensors  under  static  deformation  are: 
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(Tij  =  pdjj  -I-  -  Esjj  [2] 

where  u{  is  a  component  of  the  displacement  vector  U  = 
(u.p  u2,  u3)  in  Cartesian  coordinates  X  —  [x1}  x2,  x3),  p  is 
the  static  internal  pressure,  8^-  is  the  Kronecker  delta 
symbol,  and  E  —  E(xlf  x2;  x3)  is  the  Young’s  elastic 
modulus.  Note  in  Eq.  [l],  and  the  entire  paper,  the  lower 
index  after  a  comma  means  differentiation  with  respect 
to  the  corresponding  spatial  coordinate. 

The  static  deformation  of  the  medium  can  be  described 
by  the  equilibrium  condition: 


3 

2  +  /;  =  °.  i=  1.2.3  [3] 

7=1 

where  fj  is  the  body  force  per  unit  volume  acting  in  the  xi 
direction.  In  addition,  volume  conservation  for  an  in¬ 
compressible  medium  leads  to  the  following  relationship 
between  displacement  and  strain  components 

V  •  [/  =  El!  4"  £22  "f  e33  =  U1>1  "f  U2>2  “f  U3>3  =  0  [4] 

Using  Eqs.  [l]  and  [2]  for  stress  and  strain  components, 
and  the  incompressibility  Eq.  [4],  the  equilibrium  condi¬ 


tion  with  eliminated  internal  pressure  p  can  be  rewritten 
in  the  following  form: 

2e12(E,11  —  E,22)  4-  2(u2,2  —  um)E912  +  2e23f^i3 

—  2e13E,23  4-  (V2u2  +  co12, i)Efl  —  (V2ii!  ~  ^12*2) E,2 

4*  o)12,3Ef3  4-  V2o >i2E  4-  3(/2,i  —  fit 2)  =  0 

2e13(£,,1i  —  E,  33)  4-  2e23E,12  4-  2(u3,3  —  u1,1)Er,13 

-  2s12E,23  4-  (V2u3  +  (o13,1)E,1  4-  a)13)2E,2  [5] 

—  (V2Ui  —  &)13,3)E,3  4-  V2co13E  +  3(/3,i  —  /i,3)  =  0 


2e23(£',22  E,  33)  4-  2e13E,12  2e12E,13  +  2(u3,3 

—  u2,2)E,23  4-  (a23JlE,1  4-  (V2u3  4-  ^23  >2)  E  i2  —  (V2u2 
~  (*)23>3)E,3  +  V20)23E  +  3(/3,2  ^2*3)  =  0 

W12  =  U2>1  “  Ul>2>  ^13  =  U3>1  “  Ul>3)  W23  =  U3>2  “  U2>3 

Clearly,  the  elasticity  reconstruction  process  based  on 
Eq.  [5]  requires  accurate  measurements  of  the  displace¬ 
ment  vector,  or,  to  be  more  precise,  requires  accurate 
estimation  of  up  to  third-order  spatial  derivatives  of  the 
displacement.  Equation  [5]  can  also  be  written  in  terms  of 
spatial  derivatives  of  strain  tensor  components.  The 
unique  solution  of  the  system  of  coupled  partial  differ¬ 
ential  equations  [5]  is  determined  by  the  boundary  con¬ 
ditions,  i.e.,  the  elastic  modulus  f$(X)  must  be  given  along 
the  boundary  of  the  reconstruction  region  of  interest 
(ROI).  It  should  also  be  noted  here  that  the  analytical 
solution  of  Eq.  [5]  is  not  generally  possible,  and  numer¬ 
ical  methods  must  be  developed  to  solve  this  system  of 
partial  differential  equations. 

Based  on  the  particular  geometry  of  the  phantoms  and 
deformation  system  used  in  these  experiments,  Eq.  [5] 
was  simplified.  A  2D  approximation  of  Eq.  [5]  was  used 
because  the  imaged  plane  was  near  the  center  of  the 


FIG.  1.  Stimulated-echo  data  acquisition  and  object  deformation 
sequence.  Mechanical  transitions  occur  during  the  long  mixing 
time  (TM)  such  that  a  static  displacement  equilibrium  is  achieved. 
Local  displacement  between  deformation  states  “A”  and  “B”  are 
encoded  by  phase  shift  proportional  to  “G-displace”  amplitude 
Gd,  and  duration  t. 
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FIG.  2.  Urethane  rubber  phantom  containing 
two  hard  inclusions  and  held  between  defor¬ 
mation  plates  on  top  and  bottom  of  the  phan¬ 
tom.  Images  of  stimulated  echo  (a)  magnitude, 
and  phase  shift  proportional  to  (b)  vertical  and 
(c)  horizontal  differential  displacement.  Dis¬ 
placement  sensitivity  of  15.33  7r/mm  and  dif¬ 
ferential  displacement  of  *=*1.4  mm  between 
deformation  plates  was  used. 


C 

phantom.  For  such  a  plane  (arbitrarily  denoted  as  x3  =  0), 
the  displacement  vector  components  do  not  vary  signif¬ 
icantly  as  a  function  of  the  out-of-plane  x3  coordinate, 
and  therefore,  the  “plane  strain  state”  condition  is  appli¬ 
cable.  With  this  condition,  Eq.  [5]  reduces  to  a  single 
nontrivial  equation: 


founded  by  interference  of  the 
primary  shear  wave  with  re¬ 
flected  or  standing  shear 
waves.  To  avoid  this  condi¬ 
tion,  a  “static”  displacement 
encoding  approach  was 
adopted.  It  requires  measure¬ 
ment  of  internal  displacement 
between  two  or  multiple  de¬ 
formations  while  the  object  is 
in  mechanical  equilibrium  for 
each  measurement.  A  stimu¬ 
lated  echo  sequence  with  dis¬ 
placement-encoding  gradient 
pulses  is  used  to  achieve  this, 
as  shown  in  Fig.  1.  Mechanical 
transition  from  state  “A”  to 
state  “B”  occurs  during  the 
stimulated  echo  mixing  time, 
TM.  A  relatively  long  mixing 
time  allows  long-lived  elastic 
vibrations  to  dampen  before 
spatial  encoding.  Because  the 
relevant  magnetization  during 
TM  is  longitudinal,  it  is  unaf¬ 
fected  by  potentially  ill-de¬ 
fined  motions  during  the  me¬ 
chanical  transition  period.  As 
a  result,  a  more  accurate  static 
deformation  measurement  is 
achieved.  Also  note  that  pre¬ 
cise  synchronization  of  me¬ 
chanical  and  pulsed  gradient 
events  is  not  critical  as  long  as 
the  mechanical  transition  be¬ 
gins  after  the  second  RF  pulse 
and  is  complete  before  the  third  RF  pulse.  Similarly,  a 
long  delay  in  TE  could  be  used,  but  this  is  done  at  the 
expense  of  signal  lost  to  T2  decay. 

Local  displacement  is  encoded  by  means  of  phase  shift 
governed  by  pulsed-field  gradient  factors, 

irT^r  [7] 


(£>11  ~  £>22)(U1>2  F  U2>l)  +  2f},i2(U2>2  ~~  Ul>l) 

+  2£,,1(u2>ii  F  ^2*22)  “  2 E,2(u1,11  +  LZ i , 22)  ^ 

F  £(^2>111  F  U2>221  ~  Ul>112  ~  Ul>222)  F  3(/2>1  —  /i,2)  =  0 

Under  conditions  in  which  these  assumptions  are  valid, 
only  in-plane  displacement  u2  and  u2  or  their  derivatives 
are  needed  to  reconstruct  the  modulus  in  the  plane  x3  =  0. 

Static  Displacement  Measurement  by  Means  of 
Stimulated  Echo  MRI 

Shear-wave  propagation  speed  in  soft  tissue  is  1-20  m/s. 
Consequently,  a  shear  wave  imparted  by  a  single-stroke 
or  an  oscillating  deformation  force  may  require  tens  of 
milliseconds  to  traverse  an  object  ~100  mm  in  size.  The 
time  for  reflected  waves  to  dampen  may  be  much  longer. 
“Dynamic”  measurements,  which  encode  displacement 
during  shear-wave  propagation,  are  potentially  con¬ 


where  displacement  sensitivity,  Od,  is  in  units  of  (radi¬ 
ans/distance).  A  phase  reference  acquisition  is  required 
for  each  displacement  encode  condition  to  remove  pre¬ 
existing  phase  shifts  that  are  unrelated  to  displacement. 
Reference  data  are  acquired  using  the  same  pulse  se¬ 
quence,  including  displacement  encode  gradient  pulses, 
but  with  the  object  maintained  in  state  B.  Note  that  all 
spatial  encoding  occurs  from  the  third  RF  pulse  and 
beyond.  That  is,  the  object  is  in  one  deformation  state 
(state  B)  for  both  displacement  and  phase  reference  ac¬ 
quisitions  during  all  spatial  encoding  segments  of  the 
sequence.  Consequently,  image  registration  or  feature 
tracking  algorithms  (14)  are  not  required  to  estimate  dis¬ 
placement.  Instead,  local  displacement  is  encoded  di¬ 
rectly  by  local  phase  of  a  corrected  dataset,  Scor,  given  by, 


Scor(r) 


SA f)  SBCry 
|SB(?)| 


|SA(r)| 


[8] 
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FIG.  3.  Strain  images  calculated  from  first-order  derivatives  of  the  displacement  images  repre¬ 
sented  in  Fig.  2.  Normal  strains  (a)  e-,-,  and  (b)  e22,  and  (c)  shear  strain  s12  =  e21  reflect  internal 
elastic  properties  and  the  externally  applied  deformation  field.  The  plane  strain  state  assumption 
and  phantom  incompressibility  suggest  su  -e22  as  is  supported  by  the  relatively  featureless 
map  of  6- 1 -|  T  £22  (d). 


where  SA  and  SB  are  the  data  acquired  with  the  object 
initially  in  state  A  and  state  B,  respectively.  The  un¬ 
wrapped  phase  of  the  corrected  dataset  and  Eq.  [7]  pro¬ 
vides  a  local  measure  of  displacement,  U,  by  means  of; 

<p(?)  =  4>d  •  [r„  -  rB]  =  4>d  •  U(r )  [9] 

Most  sources  of  phase  error,  such  as  static  field  inho¬ 
mogeneity,  tend  to  be  slowly  varying  functions  of  posi¬ 
tion.  Therefore,  the  phase  reference  datasets  may  be  ac¬ 
quired  at  relatively  low  spatial  resolution  to  reduce  scan 
time. 

METHODS 

Data  Acquisition 

Elasticity  imaging  was  performed  on  two  phantoms.  One 
was  an  8 5 -mm  diameter  cylindrical  urethane  rubber 
phantom  containing  two  8-mm  cylinders  of  hard  mate¬ 
rial.  Previously,  the  ratio  of  Young’s  modulus  between 
the  inclusion  and  background  material  was  measured 


10.5  ±  1.5  for  30%  surface  de¬ 
formation  (19).  In  the  present 
study,  smaller  surface  defor¬ 
mation  was  applied,  and  con¬ 
sequently,  the  elasticity  con¬ 
trast  within  this  phantom 
should  be  less  than  the  previ¬ 
ously  measured  ratio.  Whereas 
the  mechanical  properties  of 
the  rubber  phantom  mimic 
soft  tissue,  the  phantom  had 
inherently  poor  NMR  signal.  A 
more  tissue-equivalent  phan¬ 
tom  in  terms  of  NMR  and  me¬ 
chanical  properties  was 
achieved  by  embedding  a  fresh 
canine  dog  kidney  (<24  h  ex 
vivo)  into  a  130  mm  X  105 
mm  X  75  mm  block  of  5%  gel¬ 
atin.  One  hour  before  MRI,  10 
ml  of  5%  glutaraldehyde  solu¬ 
tion  was  injected  into  kidney 
parenchyma  to  create  a  hard 
lesion(s). 

Phantoms  were  held  se¬ 
curely  in  place  under  moder¬ 
ate  preload  pressure  by  two 
parallel  acrylic  plates.  When 
pressure  to  the  top  plate  was 
released,  the  phantom  recoiled 
vertically  an  amount  con¬ 
strained  by  physical  stops. 
Maximum  vertical  displace¬ 
ment  was  <1.5  mm,  which 
represented  <2%  differential 
between  state  “A”  (greater  de¬ 
formation)  and  state  “B”  (less 
deformation).  Deformation 
was  actuated  pneumatically 
by  an  air-filled  bladder  on  top 
of  the  phantom  holder.  Pneu¬ 
matic  pressure  was  stepped  by  a  remote  solenoid  valve 
with  timing  controlled  by  an  external  transistor  transistor 
logic  (TTL)  gate  circuit  triggered  by  the  pulse  sequence. 

Displacement  encoding  gradient  pulse  duration,  r  = 
4.5  msec,  and  amplitude,  Gd  =  40  mT/m,  provided  a 
displacement  sensitivity  of  <I>d  =  15.33  7r/mm  by  means 
of  Eq.  [7].  The  displacement  encoding  direction  was  al¬ 
ternated  each  pulse  repetition  between  vertical  and  hor¬ 
izontal.  For  the  urethane  rubber  phantom,  acquisition 
parameters  were  TR  =  1.3  s,  TM  =  350  msec,  TE  =  50 
msec,  128  X  128  matrix,  four  signal  averages,  100  mm 
field  of  view,  and  10-mm  section  thickness.  An  addi¬ 
tional  128  X  32  matrix  acquisition  was  collected  for 
phase-reference  correction  of  vertical  and  horizontal  en¬ 
coded  data.  The  kidney  phantom  acquisition  parameters 
were  TR  =  1  s,  TM  =  200  msec,  TE  =  76  msec,  256  X  256 
matrix,  two  signal  averages,  150  mm  field  of  view,  and 
5-mm  section  thickness;  with  a  256  X  32  dataset  acquired 
for  phase  correction.  All  experiments  were  performed  on 
a  2  T  18-cm  bore  MRI  system  (Bruker,  formerly  GE  NMR 
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FIG.  4.  (a)  Reconstructed  map  of  Young’s  modulus  for  the  urethane 
rubber  phantom  within  a  63  X  50  mm  region.  The  boundary  of  the 
region  is  defined  to  have  a  Young’s  modulus  -  1.  Relative  Young’s 
moduli  for  inclusions  and  background  material  are  shown  in  (b). 

Instruments),  using  a  150-mm  transmit/receive  birdcage 
coil. 


Data  Processing 

Time-domain  data  were  transferred  for  off-line  process¬ 
ing  as  follows.  Phase  reference  datasets  were  zero-filled 
and  2D  Fourier  transformed  to  a  128  X  128  or  256  X  256 
matrix  for  phase  correction  by  means  of  Eq.  [7],  The 
resulting  phase  maps  were  used  to  estimate  the  spatial 
derivatives  of  the  in-plane  displacements  necessary  for 
elasticity  reconstruction  (6).  Note  phase  unwrapping  is 
not  strictly  required  because  only  phase  derivatives  are 
used.  Assuming  the  displacement  fields  are  continuous, 
resulting  in  small  differential  displacement  at  any  pixel 
compared  with  the  total  displacement,  the  differential 
displacement  between  two  neighboring  pixels  was  di¬ 
rectly  computed  from  the  angle  of  the  complex  multipli¬ 
cation  of  each  pixel  with  the  conjugate  of  the  neighboring 
pixel,  then  scaled  by  l/<l>d. 

Solving  Eq.  [6]  for  unknown  E  (xl,  x2)  performed  the 
elasticity  reconstruction,  i.e.,  reconstruction  of  the  spa¬ 
tial  distribution  of  elastic  Young’s  modulus.  As  was 
noted  previously,  the  unique  solution  of  a  boundary 
value  problem  (Eq.  [5]  or  [6])  is  determined  by  the  bound¬ 


ary  conditions.  Therefore,  a  rectangular  ROI  was  identi¬ 
fied  within  the  imaging  planes  for  both  phantoms.  For 
the  phantom  with  two  hard  inclusions,  the  ROI  was  a 
region  of  63  X  50  mm  positioned  approximately  in  the 
center  of  the  phantom  and  included  both  inclusions.  For 
the  canine  kidney  phantom,  the  rectangular  94  X  51  mm 
ROI  included  the  whole  kidney  cross-section.  In  both 
cases,  the  Young’s  modulus  value  along  the  ROI  boundary 
was  set  to  “one”  resulting  in  reconstruction  of  relative 
Young’s  modulus.  More  detailed  analysis  and  discussion  of 
defining  the  ROI  was  considered  previously  (26). 

Elasticity  reconstruction  Eqs.  [5]  and  [6]  assume  that 
spatial  derivatives  of  the  Young’s  modulus  are  continu¬ 
ous  functions.  To  ensure  continuous  elasticity  distribu¬ 
tion,  the  spatial  derivatives  of  the  displacement  were 
low-pass  filtered  before  elasticity  reconstruction,  result¬ 
ing  in  mild  spatial  resolution  reduction. 

After  defining  the  boundary  conditions,  Eq.  [6]  was 
discretized  over  the  ROI  with  the  same  grid  spacing  as 
the  MR  images,  where  all  spatial  derivatives  of  the  dis- 
placement/strain  (i.e.,  coefficients  inEq.  [6]  for  unknown 
Young’s  modulus  distribution)  were  approximated  „ by 
finite  differences.  The  linear  set  of  equations  resulting 
from  discretization  of  Eq.  [6]  was  solved  iteratively, 
where  the  error  in  each  step  was  estimated  by  averaging 
the  left-hand  side  of  Eq.  [6]  over  the  ROI  using  the  current 
estimate  of  the  elasticity  distribution.  From  step  to  step, 
the  Young’s  modulus  distribution  was  updated  based  on 
the  changes  in  the  average  error. 

RESULTS 

Magnitude  and  corrected  phase  images  of  the  urethane 
rubber  phantom  are  shown  in  Fig.  2.  Given  that  <hd  = 
15.33  7r/mm,  the  number  of  27t  phase  bands  in  Fig.  2b 
indicates  that  the  vertical  excursion  of  the  phantom  top 
relative  to  the  bottom  was  «*1.4  mm;  similarly,  the  rela¬ 
tive  lateral  displacement  between  left  and  right  edges  of 
the  phantom  was  ^0.8  mm  (Fig.  2c).  Reduced  phase 
slopes  in  the  regions  of  the  hard  inclusions  are  clearly 
visible  on  the  phase  images.  Normal  strain,  eia  and  e22, 
and  shear  strain,  e12  -  e21,  maps  are  illustrated  in  Figs. 
3a-3c,  respectively.  The  observed  contrast  reversal  be¬ 
tween  sai  and  s22  is  a  result  of  phantom  incompressibil¬ 
ity  (like  soft  tissue),  which  yields  ela  =  -e22  assuming 
negligible  out-of-plane  strain.  Consequently,  the  sum  [ett 
+  s22)  is  relatively  “flat”  as  shown  at  equivalent  gray¬ 
scale  settings  in  Fig.  3d.  Also  note,  whereas  the  strain 
maps  clearly  exhibit  object-specific  detail  (i.e.,  inclu¬ 
sions),  features  related  to  the  applied  external  deforma¬ 
tion  are  quite  conspicuous.  This  fact  demonstrates  why 
an  elasticity  reconstruction  is  needed.  The  Young’s  mod¬ 
ulus  was  reconstructed  for  a  63  X  50  mm  region  as 
presented  in  Fig.  4a.  The  boundary  of  the  elasticity  re¬ 
construction  area  was  defined  to  have  value  “1.”  The 
relative  elastic  moduli  for  select  regions  are  indicated  in 
Fig.  4b  and  are  consistent  with  the  known  elasticity  of 
these  materials  (19). 

The  canine  kidney  phantom  is  shown  in  Fig.  5  as 
magnitude  (Fig.  5a),  vertical  (Fig.  5b),  and  horizontal 
(Fig.  5c)  phase  shift  images.  It  is  apparent  from  the  ver¬ 
tical  phase  image  (Fig.  5b)  that  although  ^1.24-mm  rel- 
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FIG.  5.  Stimulated  echo  images  of  canine  kid¬ 
ney  imbedded  in  gelatin,  (a)  Magnitude  image, 
and  (b)  vertical  and  (c)  horizontal  displacement 
phase  images  in  which  <1.3  mm  of  vertical 
differential  deformation  and  15.33  7r/mm  dis¬ 
placement  sensitivity  were  applied. 
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ative  displacement  spans  the  full  phantom,  only  ^0.26 
mm  of  it  is  within  the  kidney.  Consequently,  there  is  a 
high  concentration  of  strain  near  the  gel-kidney  interface, 
as  is  clear  on  the  strain  maps  (Fig.  6).  As  noted  before, 
strain  images  exhibit  contrast  related  to  a  combination  of 
internal  structure  and  the  externally  applied  deforma¬ 
tion.  Elasticity  reconstruction,  however,  reduces  the  am¬ 
biguity  and  exhibits  contrast  dominated  by  internal  elas¬ 
tic  properties  (Fig.  7a).  It  is  also  encouraging  to  note  that 
whereas  there  was  only  moderate  strain  contrast  within 
the  kidney,  elasticity  contrast  within  renal  parenchyma 
and  central  sinus  are  well  distinguished  on  the  Young’s 
modulus  image.  Moreover,  the  site  of  glutar aldehyde  in¬ 
jection  (top-right  quadrant  of  images)  exhibited  the  high¬ 
est  relative  Young’s  modulus.  Approximately  20  h  after 
MRI,  the  kidney  phantom  was  sliced  at  a  plane  corre¬ 
sponding  to  that  studied  by  MRI.  An  optical  image  of  this 
slice  is  shown  in  Fig.  7c.  The  freshly  cut  surface  was  pal¬ 
pated  such  that  areas  of  relatively  “hard”  parenchyma 
could  be  noted.  Arrows  in  Fig.  7c  mark  the  most  conspic¬ 
uous  areas  of  hardness;  the  largest  area  corresponds  to  the 
high  Young’s  modulus  region  in  the  upper-right  quadrant 
of  the  kidney. 

DISCUSSION 

In  this  work  we  introduce  a  method  to  image  and  quan¬ 
tify  internal  elastic  properties  of  an  object  by  means  of 
displacement-sensitive  MRI  with  associated  elasticity  re¬ 
construction.  The  data  acquisition  segment  employs  gra¬ 
dient  pulses  to  encode  internal  displacement  by  means  of 


phase  using  a  stimulated  echo 
sequence.  Internal  displace¬ 
ments  occur  in  response  to  an 
external  deformation  force 
synchronized  to  the  acquisi¬ 
tion  sequence.  By  timing  con¬ 
trol,  mechanical  motion  oc¬ 
curs  while  the  relevant 
magnetization  is  longitudinal. 
The  stimulated  echo  allows 
extension  of  the  mechanical 
transition  period  to  avoid  po¬ 
tentially  long-lived  or  ill-de¬ 
fined  oscillations  within  the 
object  such  that  an  estimate  of 
“static”  displacement  is 
achieved.  Image  registration  or 
feature  tracking  (14)  is  not  re¬ 
quired  because  the  object  is  in 
one  deformation  state  for  all 
spatial  encoding.  Tt  relaxation 
and  diffusion,  which  erode 
signal,  ultimately  set  practical 
limits  on  this  period.  In  these 
experiments,  200  to  350  msec 
was  sufficient  to  allow  static 
displacement  measurement  of 
rubber  and  gelatin/tissue 
phantoms  using  a  simple  air- 
bladder  pneumatic  system. 
This  arrangement  required  the 
object  to  “passively”  recoil  during  the  TM  period. 
Clearly,  a  faster  deformation  system  can  be  built  using 
external  forces  to  “actively”  deform  the  object  during  the 
stimulated  echo  mixing  period.  A  shorter  mixing  time 
would  then  be  used  yielding  a  larger  signal. 

To  date,  two  approaches  are  present  in  elasticity  im¬ 
aging:  static  reconstructive  elasticity  imaging  (25,  26)  and 
dynamic  shear-wave  elasticity  imaging  (20,  21,  23).  In 
both,  an  external  static  or  dynamic  deformation  is  ap¬ 
plied  while  the  resulting  displacement/strain  or  propa¬ 
gating  shear  wave  is  detected  using  an  imaging  modality. 
In  reconstructive  elasticity  imaging,  the  elasticity  distri¬ 
bution  must  be  reconstructed  from  static  displacement 
and  strain  images.  The  ability  to  control  the  internal 
deformation  pattern  by  varying  the  externally  applied 
load,  and  high  SNR  displacement  and  strain  estimates 
are  the  benefits  of  this  method,  although  numerical  re¬ 
construction  algorithms  are  required.  Wherein  the  mod¬ 
eled  assumptions  are  valid,  these  algorithms  exist.  For 
more  general  applications,  they  must  be  refined  further. 
In  shear-wave  elasticity  imaging,  local  shear  wavelength 
measurements  allow  direct  and  simple  calculation  of  the 
shear  elastic  modulus.  However,  the  interference  of  shear 
waves  reflected  from  any  elasticity  inhomogeneities 
within  the  tissue,  along  with  attenuation  of  shear  waves, 
and  conversion  between  shear  and  bulk  waves  are  chal¬ 
lenges  of  this  method. 

In  these  experiments  a  displacement  sensitivity  of  <Fd 
=  15.33  7r/mm  was  achieved  using  moderate  gradient 
factors.  The  ultimate  quality  of  Young’s  modulus  recon¬ 
struction  depends  on  the  induced  phase  shift,  equal  to 
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the  product  of  <bd  and  local  displacement.  There  is  a 
limit,  however,  to  the  advantages  gained  by  increasing 
phase  shift.  As  spatial  phase  gradients  become  large,  the 
phase  distribution  within  a  given  voxel  reduces  signal 
amplitude.  Assuming  a  linear  phase  distribution  of  range 
j3  with  a  voxel,  the  signal  modulation  function  is  given  by 
sinc(/3).  For  illustration,  consider  the  vertical  phase  ex¬ 
cursion  of  21  tt  observed  across  the  80-mm  phantom  in 
Fig.  2b.  If  the  vertical  phase  excursion  was  evenly  dis¬ 
tributed  across  voxels,  the  phase  range  within  each 
0.78-mm  voxel  would  be  approximately  0.2  tt.  This  im¬ 
plies  a  signal  reduction  factor  of  sinc(0.2  ir)  =  0.94  (i.e., 
signal  loss  of  6%).  Clearly  phase  gradients  can  be  more 
concentrated  depending  on  object  geometry,  elastic  het¬ 
erogeneity,  and  deformation  geometry.  This  concentra¬ 
tion  can  lead  to  regions  of  significant  signal  loss  in  a 
manner  analogous  to  flow  dephasing  in  conventional 
MRI.  Under  such  conditions,  a  smaller  voxel  size  can 
(paradoxically)  yield  higher  signal.  In  addition,  inspec¬ 
tion  for  significant  signal  loss  within  the  displacement- 
sensitive  magnitude  image  can  identify  high-strain  re¬ 
gions  near  soft/hard  interfaces  of  a  lesion. 

Water  diffusion  in  the  presence  of  displacement  en¬ 
coding  gradients  is  another  source  of  signal  attenuation. 
We  estimate  the  signal  reduction  factor  for  freely  diffus¬ 
ing  water  was  ^0.24  in  these  experiments  (i.e.,  76% 
signal  lost  to  diffusion  effects).  Diffusion  effects  are  less¬ 
ened  by  reducing  Od,  or  alternatively  by  shortening  TM 
without  affecting  <bd.  In  either  case,  diffusion  effects  are 
assumed  independent  of  the  deformation  state,  and 


therefore  are  ignored  in  the 
elasticity  reconstruction.  Be¬ 
cause  displacement  phase 
shift  is  the  product  of  local 
displacement  and  Od,  the  se¬ 
lection  of  <i>d  is  somewhat  ar¬ 
bitrary  as  long  as  the  applied 
differential  deformation  is  ad¬ 
equate  for  elasticity  recon¬ 
struction.  In  these  preliminary 
experiments,  the  differential 
deformation  was  <1.5  mm 
across  the  imaged  object.  For 
multi-step  acquisitions,  as 
done  here,  good  reproducibil¬ 
ity  of  deformation  is  essential. 
Significant  variation  in  defor¬ 
mation  magnitude  over  the  ac¬ 
quisition  will  lead  to  phase  in¬ 
stability,  motion-like  artifact 
in  base  images,  and  errors  that 
propagate  through  the  elastic¬ 
ity  reconstruction.  It  is  a  minor 
technical  challenge  to  achieve 
relatively  high  displacement 
reproducibility  in  the  defor¬ 
mation  apparatus.  Irreproduc- 
ible  motions  that  originate 
within  the  imaged  object, 
however,  can  be  problematic 
and  are  analogous  to  undes¬ 
ired  physiologic  motion  arti¬ 
facts  in  in  vivo  diffusion  MRI.  Fortunately,  unlike  diffu¬ 
sion  and  physiologic  motions,  the  targeted  motion  in 
elasticity  imaging  is  externally  driven.  As  such,  the  dis¬ 
placement  amplitude  in  response  to  an  external  differen¬ 
tial  deformation  can  be  significantly  greater  than  irrepro- 
ducible  or  asynchronous  displacement.  For  many  in  vivo 
applications  including  the  breast,  increasing  the  differential 
deformation  severalfold  relative  to  that  applied  in  these 
phantom  studies  can  reduce  motion  artifact.  Gradient  fac¬ 
tors  and  <I>d  would  be  reduced  accordingly,  which  would 
yield  the  added  benefit  of  increased  signal  otherwise  lost  to 
diffusion  effects. 

In  practice,  the  definition  of  a  closed  contour  of  constant 
Young’s  modulus  within  the  tissue  can  be  a  challenge.  A 
hybrid  procedure  can  be  used  as  detailed  elsewhere  (26) 
and  summarized  as  follows.  The  strain  images  are  first 
processed  to  highlight  boundaries  between  regions  of  dif¬ 
ferent  elastic  modulus.  This  procedure  is  based  on  the 
stress  continuity  property  of  continuous  media  such  as 
tissue  and  can  define  regions  of  very  small  modulus  varia¬ 
tions.  After  this  boundary  detection,  closed  contours  of 
small  elasticity  variations  are  defined.  The  modulus  along 
the  contours  is  considered  constant,  thereby  providing  the 
boundary  condition  for  complete  reconstruction  of  the  elas¬ 
ticity  within  the  region  of  interest  based  on  numerical  so¬ 
lution  of  Eq.  [6].  The  elasticity  distribution  reconstructed  in 
this  way  is  the  modulus  relative  to  the  modulus  along  the 
boundary.  For  breast  elasticity  imaging,  it  is  anticipated 
that  such  a  contour  can  be  defined  a  priori  within  the 
subcutaneous  fat  that  surrounds  breast  parenchyma.  As- 
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FIG.  7.  (a)  Reconstructed  Young’s  modulus  image  within  a  rectangular  94  X  51  mm  region 
encompassing  the  kidney,  (b)  Relative  Young’s  modulus  for  the  scribed  ROIs  indicate  high  elastic 
modulus  at  the  glutaraldehyde  injection  site  in  the  upper-right  quadrant  of  the  kidney  parenchyma, 
(c)  The  optical  image  of  the  kidney  phantom  approximately  20  h  after  MRI.  Areas  marked  by  arrows 
were  noticeably  harder  as  assessed  by  sense  of  manual  touch. 


range  to  include  deep  lying  le¬ 
sions.  One  application  would 
be  measurement  of  elasticity 
in  breast  tissue  not  accessible 
to  manual  palpation.  In  situ 
studies  of  Young’s  elastic 
modulus  performed  on  sam¬ 
ples  of  breast  tissue  indicate 
that  there  is  a  large  difference 
in  elastic  modulus  between 
normal  and  pathologically 
transformed  breast  tissues. 
Others  have  analyzed  the 
Young’s  modulus  differences 
between  different  soft  tissues 
and  have  found  1-2  orders  of 
magnitude  difference  in 
Young’s  elastic  moduli  of  a  tis¬ 
sue  in  different  physiologic 
states  (1).  If  elastic  changes 
predate  formation  of  calcifica¬ 
tions,  elasticity  imaging  could 
potentially  increase  detection 
and/or  characterization  of  ma¬ 
lignant  breast  masses  and  thus 
be  an  important  addition  to  ex¬ 
isting  clinical  diagnostic  tools. 
Practical  issues  such  as  the  rel¬ 
atively  high  cost  of  MRI  may 
hinder  use  of  this  approach  as  a 
screening  test.  Nevertheless,  ad¬ 
ditional  work  to  define  the  role 
of  this  technique  as  a  primary 
diagnostic  tool  or  supplemental 
problem-solving  modality  in 
the  management  of  soft-tissue 
disease  is  well  justified. 


suming  the  Young’s  modulus  of  the  fat  boundary  is  con¬ 
stant,  the  relative  Young’s  modulus  image  of  parenchyma 
can  be  reconstructed.  Alternatively,  the  breast  can  be  sur¬ 
rounded  by  a  high-signal  cuff  of  known  elastic  modulus 
and  imaged.  An  image  of  absolute  Young’s  modulus  can 
then  be  reconstructed  if  the  boundary  contour  is  defined 
within  the  cuff  material. 

Some  artifacts  present  in  elasticity  images  (Figs.  4  and  7) 
are  due  to  violation  of  the  plane  strain  state  approximation 
in  these  experiments.  Indeed,  if  a  plane  strain  state  is  not 
present,  the  reconstruction  based  on  Eq.  [6]  will  be  in  error. 
The  elasticity  reconstruction,  however,  does  not  have  to  be 
limited  by  a  plane  strain  assumption  if  all  3D  components 
of  the  displacement  vector  are  available.  Fortunately,  such 
information  would  be  available  using  3D  displacement  en¬ 
coding  within  a  volumetric  imaging  sequence.  The  issue  of 
long  scan  time  could  be  resolved  by  incorporating  echo- 
planar  imaging  or  fast-spin-echo  segments  for  spatial  en¬ 
coding.  Correspondingly,  an  elasticity  reconstruction  based 
on  Eq.  [5]  would  be  applied  to  produce  volumetric  elastic¬ 
ity  maps. 

The  long-range  goal  of  quantitative  elasticity  imaging 
is  to  provide  remote  palpation  thus  expanding  its  limited 
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introduction 

Elasticity  MRI  is  the  reconstruction  of  the  elastic  modulus  in 
media  using  local  displacement  measurements  of  an  object  under 
static  deformation  [1],  or  directly  from  measurements  of  shear 
wave  propagation  [2].  Often  it  is  assumed  out  of  plane  strain  is 
approxiamtely  zero.  However,  this  plane  strain  state 
approximation  does  not  hold  in  general.  Estimates  of  the  full  3D 
displacement  and  strain  fields  are  required  to  accurately 
reconstruct  elasticity.  Here  we  describe  such  a  method  based  on 
static  displacement,  stimulated-echo  imaging  [3].  As  previously 
described,  this  approach  encodes  local  displacement  in  response 
to  an  externally-applied  deformation.  A  relatively  long  STE 
mixing  time  (200-300msec)  allows  ill-defined  mechanical 
vibrations  to  dampen  to  where  equilibrium  conditions  apply.  The 
approach  uses  pulsed-field  gradients  (PFG)  to  encode 
displacement  and  is,  in  principle,  readily  extended  to  measure  the 
full  3D  displacement  field. 

METHODS 

Features  of  the  displacement-encoding  STE  technique  have  been 
described  previously  [3]  but  are  summarized  here.  The  initial 
spatial  configuration  of  an  elastic  object  is  encoded  by 
application  of  a  PFG  between  the  first  two  non-selective  90°  rf 
pulses.  An  externally-applied  deformation  force  transforms  the 
object  to  a  new  configuration  during  the  mixing  time,  TM,  of  the 
STE.  The  differential  surface  deformation  can  be  relatively 
subtle  (<5mm)  and  is  chosen  to  complement  the  applied  PFG 
area  which  defines  displacement  sensitivity  (~3-167t/mm). 
Mixing  time  must  be  sufficient  for  the  object  to  come  to  rest  by 
the  3rd  slice/slab-selective  90°  rf  pulse.  A  second  PFG  prior  to 
signal  readout  yields  a  phase  directly  proportional  to  local 
displacement  between  initial  and  final  object  configuration.  To 
achieve  3D,  a  fast-spin-echo  train  was  appended  to  the  STE  with 
an  additional  centric-ordered  z-slab  phase-encode  echo  train  (8- 
16  ETL  on  kz).  Conventional  phase-encoding  was  used  for  the 
remaining  spatial  dimension.  Acquisition  of  data  sets  with 
inteleaved  PFG  gradients  applied  along  X,  Y,  and  Z  axes: 
TR=T';  256  x  128  x  8;  4ave;  3.8Ji/mm  sensitivity.  Phase 
correction  was  done  via  a  sparse  reference  dataset  (32  ky  lines). 


A  tissue-mimicking  block  phantom  containing  ramped  bars  of 
harder  material  (~  6x  Young’s  modulus)  was  held  between  two 
pneumatic  deformation  plates.  These  parallel  plates  released  a 
5mm  differential  surface  deformation  during  TM=270ms. 

RESULTS 

Images  of  the  phantom  shown  in  figure  1  are  (a)  magnitude  and 
phase  (i.e.  displacement)  images  along  (b)  X,  (c)  Y  and  (d)  Z 
directions.  Corresponding  strain  images  illustrated  in  figure  2 
are  (a)  £xx  ,  (b)  Eyy  ,  (c)  txy,  and  (d)  (Exx+Eyy)*  Note,  (Exx+Eyy) 
should  be  zero  if  the  plane-strain  condition  holds.  Apparent 
structure  in  figure  2(d)  indicates  there  is  out  of  plane  strain. 

DISCUSSION 

Extension  of  elasticity  MRI  to  3D  is  essential  to  implement  a 
general  elasticity  reconstruction.  3D  elasticity  reconstruction 
algorithms  have  been  described  [1],  but  have  not  been 
implemented  due  to  lack  of  suitable  3D  data  acquisition  schemes. 
The  method  presented  here  represents  an  initial  step  toward  that 
end. 


(c)  £xy  (d)  6xx"b  Eyy 

Figure  2 
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Abstract.  This  article  presents  a  method  for  measuring  three-dimensional 
mechanical  displacement  and  strain  fields  using  stimulated  echo  MRI.  Additional 
gradient  pulses  encode  internal  displacements  in  response  to  an  externally  applied 
deformation.  By  limiting  the  mechanical  transition  to  the  stimulated  echo  mixing 
time,  a  more  accurate  static  displacement  measurement  is  obtained.  A  three- 
dimensional  elasticity  reconstruction  within  a  region  of  interest  having  a  uniform  shear 
modulus  along  its  boundary  is  performed  by  numerically  solving  discretized  elasticity 
equilibrium  equations.  Data  acquisition  and  reconstruction  were  performed  using  a 
silicone  gel  phantom  containing  an  inclusion  of  known  elastic  properties.  A  comparison 
between  two-dimensional  and  three-dimensional  reconstructions  from  simulated  and 
experimental  displacement  data  shows  higher  accuracy  from  the  three-dimensional 
reconstruction.  The  long  term  objective  of  this  work  is  to  provide  a  method  for 
remotely  palpating  and  elastically  quantitating  manually  inaccessible  tissues. 
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1.1.  Motivation 

Palpation  has  long  been  used  by  physicians  as  a  means  to  detect  disease.  The  underlying 
basis  for  this  detection  is  the  presence  of  “hard”  tissue.  Evidence  suggests  that  Young’s 
(or  shear)  elastic  moduli  may  differ  by  orders  of  magnitude  within  soft  tissues  in  various 
physiologic  states  (Sarvazyan  et  al  1995),  (Skovoroda  et  al  1995b).  In  addition,  manual 
self-examination  is  the  first  diagnostic  line  of  defense  against  both  breast  (Hill  et  al 
1988),  (Newcomb  et  al  1991)  and  testicular  cancers.  With  breast  cancer,  manual 
detection  of  a  new  mass  often  merits  excisional  biopsy,  even  if  uncorroborated  by  other 
tests,  as  nodule  hardness  raises  suspicion  of  malignancy  (Foster  1996).  Palpation  of 
superficial  lymph  nodes  and  abdominal  organs  is  also  routinely  performed.  Although  the 
touch  of  a  skilled  physician  is  a  powerful  diagnostic  tool,  palpation  sensitivity  is  relatively 
poor  within  deep,  dense  or  heterogeneous  tissue.  Thus,  most  manually  detected  lesions 
are  either  superficial,  relatively  large,  or  both. 

1.2.  Elasticity  imaging 

Currently,  many  scientists  are  working  on  extending  the  range  and  sensitivity  of 
palpation  by  using  various  methods  to  image  tissue  elasticity.  The  basic  method  for 
creating  an  elasticity  map  involves  two  steps.  First,  the  internal  displacements  within 
tissue  under  an  applied  mechanical  stress  are  measured.  The  (usually  externally)  applied 
deformation  may  be  either  dynamic  or  static.  Then,  from  these  data,  a  reconstruction 
of  regional  variations  in  tissue  elasticity  is  performed,  either  directly  or  after  calculating 
internal  strains.  Although  both  internal  displacements  and  strains  are  related  to  the 
elastic  properties  of  tissue,  both  are  strongly  affected  by  geometry.  Thus,  some  form  of 
reconstruction  is  necessary  to  uniquely  determine  the  elasticity  distribution. 

To  date,  two  major  medical  imaging  modalities  have  been  used  to  measure  tissue 
displacement:  ultrasound  and  magnetic  resonance  imaging  (MRI).  The  phase  sensitivity 
of  these  methods  lends  itself  to  tracking  tissue  motion.  Most  elasticity  imaging  has  been 
carried  out  using  ultrasonically  measured  tissue  displacements.  These  data  have  been 
obtained  by  tracking  specular  reflections  (Dickinson  and  Hill  1982),  (Tristam  et  al  1986, 
1988),  by  Doppler  techniques  (Lerner  et  al  1990),  (Parker  et  al  1990),  (Parker  and  Lerner 
1992),  by  cross-correlation  of  acoustic  echoes  (Ophir  et  al  1991),  (Garra  et  al  1997),  and 
by  speckle  tracking  (Adler  et  al  1988),  (O’Donnell  et  al  1994),  (Emelianov  et  al  1995a). 
Other  efforts  employ  MRI  for  measuring  tissue  motion,  as  discussed  below. 

1.3.  MRI  measurement  of  tissue  displacement 

In  the  past,  myocardial  motion  and  strain  have  been  measured  using  spatial 
magnetization  tagging  (Axel  and  Dougherty  1988),  (Zerhouni  et  al  1988),  and  phase- 
based  velocity  encoding  (Pelc  et  al  1995).  More  recently,  methods  have  been  devised 
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to  measure  tissue  displacement  specifically  for  elasticity  imaging.  These  measurements 
can  be  separated  based  upon  the  nature  of  the  applied  deformation. 

1.3.1.  Dynamic  deformation  With  these  methods,  a  periodic  excitation  is  applied  to 
the  tissue  near  the  region  of  interest,  and  the  entire  system  may  be  allowed  to  reach 
steady  state.  One  or  several  “snap  shots”  of  mechanical  wave  propagation  within  the 
object  are  produced  by  controlling  the  relative  phase  between  the  mechanical  excitation 
and  the  motion-encoding  gradients.  The  local  displacement  information  in  these  images 
is  then  used  as  an  input  for  an  elasticity  reconstruction  algorithm.  Initial  experiments 
used  a  shear  excitation,  and  the  elasticity  reconstruction  was  performed  assuming  the 
recorded  image  contained  only  shear  waves  (Muthupillai  et  al  1995).  If  only  shear  waves 
are  present  in  a  purely  elastic  medium,  local  elastic  modulus  variations  are  determined 
via  the  relation  p  =  u2X2 p,  where  p  is  the  local  shear  modulus,  v  is  the  frequency 
of  the  applied  deformation,  A  is  the  measured  local  strain-wave  wavelength,  and  p 
is  the  density  of  the  medium.  Although  attractive  in  its  simplicity,  this  approach  is 
compromised  by  frequency  dependent  visco-elastic  effects  and  strain-wave  wavelength, 
interference  from  reflections  off  of  elastic  inhomogeneities,  and  the  possible  presence  of 
longitudinal  mechanical  waves  in  the  medium.  Despite  these  limitations,  this  method 
has  been  applied  in  vivo  (Dresner  et  al  1999),  (Lawrence  et  al  1999).  Recently,  a  more 
general  elasticity  reconstruction  from  a  series  of  “instantaneous”  steady  state  mechanical 
wave  images  has  been  developed  (Sinkus  et  al  1999).  This  and  another  technique  (Van 
Houten  et  al  1999),  (Weaver  et  al  1999)  rely  on  a  more  complete  visco-elastic  tissue 
model  than  that  presented  in  Muthupillai  et  al  (1995). 

1.3.2.  Static  deformation  Another  method  of  producing  an  internal  strain  field  in  an 
object  is  to  deform  it  and  allow  the  material  to  relax  to  equilibrium  before  measuring  the 
displacement  field.  The  displacement  field  has  been  accessed  using  spatial  magnetization 
tagging,  but  this  method  suffers  from  spatial  resolution  limited  by  the  tagged  grid  size 
and  typically  measures  only  two-dimensional  (2D)  motion  (Fowlkes  et  al  1995).  A 
quasi-static  method  using  bipolar  gradient  phase  encoding  of  2D  motion  is  presented 
by  Plewes  et  al  (1995,  1996).  Stimulated  echo  MRI  has  also  been  used  to  measure  2D 
displacement  fields,  from  which  elasticity  images  have  been  reconstructed  (Chenevert 
et  al  1998).  This  method  has  been  extended  to  study  myocardial  motion  (Aletras  et 
al  1999).  With  these  techniques,  visco-elastic  effects  are  generally  ignored,  making  the 
reconstruction  more  straight-forward.  Care  must  be  taken,  however,  to  justify  the  use  of 
a  static  model,  especially  when  repeated  deformations  are  needed  to  acquire  a  complete 
data  set. 

In  general,  MRI  has  several  advantages  over  ultrasound  with  respect  to  elasticity 
imaging.  Although  ultrasound  accurately  measures  motion  along  the  beam  axis,  lateral 
motion  is  measured  with  a  resolution  given  by  the  depth-dependent  beam  width.  Out- 
of-plane  motion  is  generally  not  considered  given  the  problems  with  three-dimensional 
(3D)  image  registration  in  ultrasound.  These  restrictions  compromise  the  quality  of 
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displacement  data  available  and  constrain  the  type  of  model  used  to  produce  an  elasticity 
image.  Ultrasound  does,  though,  offer  the  advantages  of  low  cost  and  real-time  imaging. 
MRI,  on  the  other  hand,  gives  one  the  ability  to  measure  3D  displacements  within  an 
object,  and  does  this  at  a  higher  overall  resolution  than  clinical  ultrasound. 

In  this  paper  we  present  a  method  for  encoding  the  full  3D  displacement  field  within 
an  object  undergoing  an  externally  applied  static  (or  quasi-static)  deformation.  Local 
strain  estimates  are  calculated  from  the  measured  displacements,  and  the  strain  tensor 
is  used  to  numerically  solve  differential  elasticity  equilibrium  equations,  ultimately 
producing  a  3D  elasticity  image. 

2.  Reconstructive  elasticity  imaging  from  static  displacement  fields 

The  goal  of  elasticity  imaging  is  to  produce  a  map  of  the  tissue  elastic  modulus  in  a  region 
of  interest  using  available  measurements  of  displacement  components.  In  this  work,  the 
reconstruction  approach  taken  is  based  upon  a  model  of  linear,  elastic,  isotropic  media 
(Emelianov  et  al  1995b),  (Skovoroda  et  al  1995a,  1999).  The  central  equations  and 
concepts  are  covered  briefly  here.  A  more  detailed  discussion  can  be  found  in  the 
references  mentioned. 


2.1.  Linear  elasticity  and  reconstruction 


In  linear  elasticity,  the  components  of  the  strain  (e8j)  and  stress  (cr,j)  tensors  in  an 
incompressible  medium  undergoing  small  deformations  are  given  by: 


&ij  —  pSij  +  2  p£{j  (2) 


where  u,  is  a  component  of  the  displacement  vector  U  =  (uj,  it  2,  u3)  in  Cartesian 
coordinates  r  =  (xi,X2,X3),  p  is  the  product  AV  •  U  for  compressible  media  or  the 
static  internal  pressure  for  incompressible  media,  5ij  is  the  Kronecker  delta  function,  A 
and  p  are  the  Lame  coefficients,  and  p  =  p(r)  is  the  shear  elastic  modulus. 

A  medium  undergoing  static  deformation  obeys  the  equilibrium  condition: 

+  <  =  1’2’3’  (3) 

j=i  J 


where  /,•  is  the  body  force  per  unit  volume  acting  in  the  x ;  direction.  In  addition,  if  a 
medium  is  incompressible,  volume  conservation  leads  to  the  following  relation: 


,r  dui  du2  du3 

U  —  £n  +  £22  +  £33  —  7J - b  T. - b  T. —  —  0. 

OX  i  OX  2  OX  3 


(4) 


Although  not  necessary  in  the  development  that  follows,  it  should  be  noted  that  soft 
tissue  is  approximately  incompressible. 

Using  equations  (1)  and  (2)  in  (3),  the  unknown  p(r)  can  be  eliminated  to  yield  a  set 
of  differential  equations  depending  only  on  U ,  first  and  higher-order  spatial  derivatives 
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of  U ,  and  the  elasticity  distribution,  p{r).  This  set  of  equations  is  then  numerically 
solved  to  estimate  the  unknown  shear  elasticity  distribution. 


2.2.  Importance  of  three-dimensional  reconstruction  methods 


Several  approaches  have  been  proposed  to  estimate  tissue  elasticity  from  the 
experimentally  measured  spatial  distribution  of  internal  displacements  within  an  object. 
The  simplest  method  is  a  one-dimensional  (ID)  estimation  of  normalized  tissue  elasticity, 
expressed  as: 


«i  =  1/e,  (5) 

where  e  is  longitudinal  strain  (Ophir  et  al  1991),  (Garra  et  al  1997).  Indeed,  a  loaded 
object  generally  exhibits  low  longitudinal  strain  in  relatively  hard  regions  and  high 
longitudinal  strain  in  relatively  soft  regions. 

A  2D  elasticity  reconstruction,  based  on  a  plane-strain  assumption  and  all  necessary 
in-plane  strain  components,  provides  a  more  accurate  representation  of  the  object’s 
elasticity  (Skovoroda  et  al  1995a,  1999).  The  theory  of  reconstructing  clearly  bounded 
and  spatially  distributed  tissue  inhomogeneities  has  been  demonstrated  by  Skovoroda 
et  al  (1995a)  as  well.  However,  inaccurate  estimates  may  result  by  using  either  a  ID  or 
2D  reconstruction  of  a  3D  object. 

To  demonstrate  these  inaccuracies,  consider  a  spherical  inclusion  of  radius  R  in 
a  uniaxially,  uniformly  loaded,  infinite,  homogeneous  medium  (Goudier  1933).  For 
an  incompressible  medium,  the  distribution  of  longitudinal  strain  along  the  x3  axis 
(orthogonal  to  the  applied  deformation),  is  (Skovoroda  et  al  1994): 

■ To  <  /? 

3+2ko  ^3  —  n 

e=<  /?{1  +  «S&y[5(ft)3-9(S)5]}  "  (6) 


x3  >  R. 


Here  j3  is  the  magnitude  of  the  applied  strain  and  «o  =  p/po  is  the  ratio  of  the  inclusion 
to  background  shear  moduli. 

Normalizing  (6)  by  /?,  which  corresponds  to  the  axial  strain  in  the  tissue  far  from 
the  inclusion,  and  substituting  into  (5)  we  obtain: 


3+2kq 

5 


x3  <  R 


« 1  = 


(3  +  2/c0)  ^3  -f  2k0  -f  5  (^-)  -9^^  j  x3>R. 


(7) 


Note  that  /ci/ko  =  (3  +  2k0)/5/co  within  the  inclusion.  That  is,  for  a  very  hard  inclusion 
(k0  large),  the  relative  modulus  obtained  from  a  ID  reconstruction  will  only  be  40% 
of  its  actual  value.  On  the  other  hand,  for  a  soft  inclusion  (k0  small),  the  relative 
modulus  estimate  will  approach  3/5,  no  matter  how  much  softer  the  inclusion  is  than 
the  background.  Obviously,  the  inaccuracy  of  a  ID  elasticity  estimation  may  not  be 
acceptable  for  many  applications. 

Now  consider  a  2D  reconstruction.  Figures  1(a)  and  (e)  show  the  exact  relative 
elasticity  distribution,  k(x i,ar2),  for  two  infinitesimal  planes  in  our  pedagogic  phantom 
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with  kq  =  4.  Figure  1(a)  presents  k  for  x3/R  ~  1.05,  that  is,  outside  of  the 
inclusion,  while  figure  1(e)  is  the  x3/R  sa  0.95  plane.  The  corresponding  relative  2D 
reconstructions,  /c2(a?i,  X2),  are  shown  in  figures  1(b)  and  (f).  The  reconstructions  were 
performed  using  the  algorithm  presented  by  Skovoroda  et  al  (1999).  For  comparison 
with  experimental  results  (see  section  5),  an  analytic  model  was  used  to  generate 
displacement  data  which  was  sampled  with  the  x2  resolution  of  the  experimental 
displacement  encoded  data  discussed  in  section  4.2.  The  strains  used  as  input  for  the 
reconstructions  were  calculated  as  described  in  section  4.3,  and  the  reconstructions  were 
performed  over  a  region  of  interest  identical  to  the  one  discussed  in  that  same  section. 
The  positions  of  the  two  reconstructed  planes  were  selected  to  apporximately  correspond 
to  the  experimental  planes  presented  in  section  5.  As  evidenced  here,  neglecting  out- 
of-plane  strain  components  in  the  reconstruction  results  in  geometrical  distortions  in 
the  elasticity  image.  In  this  case,  the  spherical  inclusion  is  reconstructed  as  a  prolate 
spheroid.  The  reconstruction  inaccuracy  of  a  plane-strain  based  reconstruction  is  small 
near  the  central  plane  and  increases  with  the  distance  between  the  imaging  plane  and 
the  center  of  the  inclusion.  Far  from  the  inclusion,  a  2D  reconstruction  would  again  be 
accurate. 

It  is  clear  that  a  ID  or  2D  reconstruction  may  lead  to  significant  inaccuracies  in 
tissue  elasticity  estimations,  especially  when  complicated  in  vivo  geometries  influence 
displacement  and  strain  measurements.  This  points  to  the  need  for  an  accurate  3D 
elasticity  reconstruction.  As  shown  in  (Skovoroda  et  al  1995a,  1999),  a  general  unknown 
shear  elasticity  distribution,  fi(x  1,0:2,  £3),  must  satisfy  the  equation: 

dVen)  _  d2Quei2)  d2[/*(e2 2  -  en)]  gVfga)  dfygis)  _  ,  , 

dx\  dx\  dx\dx2  dx\dx3  dx2dx3 

Thus,  in  order  to  compute  all  the  necessary  components  of  the  strain  tensor,  e,j,  in 
(8),  all  of  the  displacement  components  (ui,u2,  u3)  must  be  measured  as  a  function  of 
spatial  coordinates  (aq,  x2,  x3).  This  requirement  exists  in  both  the  differential-based  3D 
reconstruction  (8),  as  well  as  in  the  more  stable  integral  based  3D  approach  (Skovoroda 
et  al  1999). 

The  3D  elasticity  reconstructions  from  the  two  planes  previously  discussed  are 
shown  in  figures  1(c)  and  (g).  The  reconstruction  was  performed  as  discussed  in 
section  4.3.  Although  not  perfect  due  to  the  relatively  large  x3  step  size,  the  3D 
reconstructions  clearly  exhibit  fewer  geometric  distortions  than  the  2D  estimates.  This 
is  particularly  well  illustrated  by  the  central  vertical  profiles  through  the  analytic,  2D, 
and  3D  shear  distributions  presented  in  figures  1(d)  and  (h).  In  the  x3/R  sa  1.05  plane, 
the  2D  reconstruction  estimates  that  an  inclusion  is  present,  when  indeed  it  is  not, 
while  the  3D  reconstruction  shows  little  evidence  of  the  presence  of  an  inclusion.  The 
estimate  of  the  extent  of  the  inclusion  in  the  x3/R  sa  0.95  plane  is  also  improved  over 
the  2D  estimate.  As  with  the  2D  reconstructions,  the  strain  data  and  reconstruction 
parameters  used  for  the  3D  reconstruction  were  identical  to  those  of  the  experimental 
parameters  described  in  sections  4.2,  4.3,  and  5. 
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(a)K(xrx2)  (b)K2(xrx2) 


(C)  K3(XrX2) 


(d)  Central  profiles 
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Figure  1.  Simulated  elasticity  distributions,  k(xi,  X'2),  and  corresponding  2D, 
ki(xi,X2),  and  3D,  K3(®i,  x2)  elasticity  reconstructions  from  the  x$ /R  ss  1.05,  (a)-(d), 
and  xz/R  ~  0.95,  (e)-(h),  planes  of  a  phantom  with  a  single  hard,  spherical  inclusion 
of  radius  R.  Also  presented  are  the  central  vertical  profiles  of  each  distribution,  where 

( - )  is  k,  ( . •)  is  /C2,  and  ( - )  is  k-j.  All  are  presented  on  a  log  scale  where 

black  corresponds  to  a  relative  shear  modulus  of  0.5  and  white  to  4.5.  The  background 
has  a  realtive  shear  modulus  of  1,  and  the  inclusion,  4. 


3.  Static  displacement  measurement  via  stimulated  echo  MRI 

As  previously  discussed,  static  displacement  measurements  for  elasticity  imaging  avoid 
several  confounding  factors  that  may  be  present  if  dynamic  displacement  measurements 
are  used.  Since  shear  wave  propagation  speed  in  soft  tissue  is  approximately  1-20  m  s-1, 
shear  waves  launched  into  a  medium  by  an  applied  deformation  may  require  tens  of 
milliseconds  to  traverse  an  object  approximately  100  mm  in  size.  Reflected  waves 
may  take  much  longer  to  dampen.  To  appropriately  measure  an  object’s  internal 
static  displacements,  the  object  must  be  in  mechanical  equilibrium — that  is,  it  must 
satisfy  (3) — during  both  the  pre-  and  post-deformation  measurements.  A  stimulated 
echo  MRI  sequence  using  displacement  encoding  gradient  pulses  is  employed  to  achieve 
this  (Chenevert  et  al  1998).  Figure  2  presents  a  schematic  of  this  pulse  sequence. 
The  mechanical  transition  from  the  pre-  to  post-deformational  states  occurs  during  the 
stimulated  echo  mixing  time,  Tm .  Because  the  relevant  magnetization  is  longitudinal 
during  Tm ,  it  is  unaffected  by  the  object’s  motion  during  the  mechanical  transition 
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Figure  2.  Displacement  encoding,  stimulated  echo  pulse  sequence  waveforms.  RF  = 
radio  frequency,  Ga  =  displacement  encoding  gradient,  and  Gro  =  read-out  (a?i),  Gpe  = 
phase-encode  (£2),  and  Gsi  =  slice  (2:3)  directed  gradient  waveforms.  Tm  is  the  mixing 
time,  Te  is  the  echo  time,  and  r  is  the  duration  of  the  displacement  encoding  gradient. 
Note  that  the  displacement  encoding  gradient  may  be  applied  to  any  of  the  directional 
waveforms. 


period.  This  allows  a  more  accurate  measurement  of  static  internal  displacement. 
Additionally,  precise  synchronization  of  the  motion  and  applied  gradients  is  not 
necessary  as  long  as  the  mechanical  deformation  begins  after  the  second  radio-frequency 
pulse,  and  internal  motion  stops  before  the  third.  A  long  delay  in  the  echo  time,  T#, 
could  also  be  used  to  let  the  object  reach  equilibrium,  but  this  would  likely  lead  to 
prohibitive  signal  loss  from  T2  decay. 

Local  displacements  are  encoded  in  the  magnetization’s  phase  via  pulsed-field 
gradients.  The  displacement  sensitivity,  in  radians/distance,  of  the  sequence  is: 


=  7  [  Gd(t)dt  -  7 Gdr, 
Jo 


(9) 


where  7  is  the  gyromagnetic  ratio  of  the  proton,  Gd(t)  is  the  encoding  gradient 
waveform,  and  r  is  the  duration  of  the  encoding  gradient.  However,  for  accurate 
displacement  measurements,  phase  shifts  unrelated  to  the  applied  deformation  must 
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be  removed.  This  is  done  by  acquiring  a  phase  reference  data  set  using  the  same  pulse 
sequence,  but  with  the  object  maintained  in  the  post-deformational  state  for  the  entire 
experiment.  Note  that  all  spatial  encoding  takes  place  with  the  object  in  the  post- 
deformational  state  for  both  the  displacement  encoded  and  reference  acquisitions.  Thus, 
no  image  registration  or  tracking  algorithms  are  required  to  use  the  reference  data,  Sr, 
to  correct  the  displacement  encoded  data,  S<i.  The  corrected  data  set,  Sc,  is  then: 

W  =  *  l&M|e*|r)-  (10) 

Most  sources  of  phase  error,  such  as  static  field  inhomogeneities,  tend  to  be  slowly 
varying  functions  of  position.  Thus  the  phase  reference  data  may  be  acquired  at 
relatively  low  spatial  resolution  to  reduce  scan  time. 

The  unwrapped  phase  of  (10)  is  related  to  the  local  displacement  vector,  U ,  via: 

4>(r)  =  <&d  •  Ar  =  3>j  •  U(r),  (11) 

where  A r  is  the  local  displacement  from  pre-  to  post-deformational  states.  The 
displacement  sensitivity,  may  be  made  sensitive  to  motion  in  an  arbitrary  direction 
based  upon  appropriate  combination  of  displacement  encoding  gradients  in  the  read¬ 
out,  phase-encode,  and  slice  directions.  Hence,  this  pulse  sequence  readily  extends  to 
acquiring  three-dimensional  displacement  data. 

4.  Methods 

4-1.  Phantom 

Elasticity  imaging  experiments  were  performed  on  a  phantom  with  a  spherical  hard 
inclusion.  Semicosil  921  silicone  gel  (Wacker  Silicones  Corporation,  Adrian,  MI)  was 
used  to  construct  a  phantom  qualitatively  simulating  the  mechanical  properties  of  soft 
tissue.  The  Semicosil  921  consists  of  two  components,  A  and  B,  wherein  different  ratios 
of  these  components  are  used  to  vary  the  mechanical  properties  of  the  gel.  A  tissue- 
mimicking  phantom  was  constructed  in  several  steps.  First,  background  material  was 
prepared  by  thoroughly  mixing  components  A  and  B  in  a  1:1  ratio,  and  then  pouring 
the  mixture  into  a  154-mm  by  80-mm  rectangular  mold.  The  mixture  was  degassed  and 
cured  for  24  hours  at  room  temperature  to  produce  a  22-mm  thick  layer.  Then  a  25-mm 
diameter  hard  sphere  was  prepared  from  a  1:2.5  mixture  of  A  and  B  and  was  placed 
on  top  of  the  layer  in  the  middle  of  the  mold.  Finally,  another  batch  of  background 
material  (1:1  ratio)  was  poured  into  the  mold  resulting  in  a  64-mm  by  80-mm  by  154- 
mm  phantom  with  a  single,  hard,  spherical  inclusion  roughly  in  the  center.  At  the 
same  time,  three  samples  of  each  batch  were  taken  to  independently  assess  the  elasticity 
contrast  between  the  inclusion  and  surrounding  materials.  These  measurements  were 
performed  using  the  force-deformation  system  described  in  (Erkamp  et  al  1998),  and 
showed  that  the  inclusion  was  four  times  harder  than  the  background,  and  that  both 
background  materials  were  elastically  equivalent. 
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4-2.  Data  acquisition 

To  provide  repeatable  deformation,  the  phantom  was  placed  under  moderate  pre-load 
pressure  between  two  acrylic  plates  in  a  pneumatically  driven  device.  Air-filled  neoprene 
boots  in  a  push-push  configuration  provided  the  necessary  force  to  the  top  plate  to  keep 
the  phantom  in  this  pre-load  state,  and  aided  the  vertical  recoil  of  the  phantom  to  the 
post-deformation  state.  Pneumatic  pressure  was  delivered  via  two  solenoid  valves  whose 
timing  was  controlled  by  an  external  transistor-transistor  logic  circuit  triggered  by  the 
pulse  sequence.  Quick-release  valves  aided  in  depressurizing  the  boots.  Both  the  pre¬ 
load  and  recoil  positions  of  the  top  acrylic  plate  were  set  by  adjustable  stops;  the  bottom 
plate’s  position  was  fixed.  The  applied  vertical  deformation  was  approximately  2.4  mm, 
or  about  6%  strain,  between  the  pre-transition  (greater  deformation)  and  post-transition 
(less  deformation)  states. 

During  data  acquisition,  the  displacement  encoding  gradient  pulse  duration,  r, 
was  1.5  ms,  and  the  amplitude,  G<\,  was  40  mT  m-1  in  the  read-out  (®i)  and  phase- 
encode  (#2)  directions,  and  60  mT  m-1  in  the  slice  (£3)  direction.  Here,  the  £3 
direction  was  along  the  bore’s  axis,  and  the  £1  and  £2  directions  were  perpendicular 
to  £3  in  the  horizontal  and  vertical  directions,  respectively.  By  (9),  the  displacement 
sensitivity,  was  approximately  5.11  ?r  mm-1  in  the  x\  and  £2  directions,  and  about 
7.66  7r  mm-1  in  the  £3  direction.  The  displacement  encoding  direction  was  cycled  each 
pulse  repetition  between  the  £1,  £2,  and  £3  directions.  The  pulse- to-pulse  repetition 
time  was  approximately  0.98  s,  the  mixing  time  ( Tm )  was  270  ms,  and  the  echo  time 
(Te)  was  45  ms.  Two  averages  were  taken  of  a  256  x  256  x  32  matrix  covering  an  80-mm 
by  110-mm  by  48-mm  field  of  view.  The  phase  reference  data  were  collected  using  a 
256  x  32  x  32  matrix  while  keeping  all  other  parameters  the  same.  All  experiments  were 
performed  on  a  2  T,  18-cm  bore  MRI  system  (Bruker,  formerly  GE  NMR  Instruments) 
using  a  150-mm  transmit/receive  birdcage  coil. 

4-3.  Data  processing 

All  time-domain  data  were  transferred  off-line  for  processing.  For  phase  correction,  the 
phase  reference  data  set  was  zero-filled  to  a  256  x  256  x  32  matrix.  Then  this  and 
the  displacement  encoded  data  were  3D  Fourier  transformed  and  corrected  as  in  (10). 
The  resulting  phase  maps  were  then  used  to  estimate  the  spatial  derivatives  to  compute 
the  strains,  via  (1),  necessary  for  the  elasticity  reconstruction.  Phase  unwrapping  of 
the  displacement  data  was  not  required  since  only  phase  derivatives  were  used  in  the 
strain  calculations.  The  displacement  derivative  at  the  ith  point  in  the  j  direction  was 
computed  from  the  angle  of  the  complex  multiplication  of  the  i  +  1th  point  with  the 
conjugate  of  the  i  —  1th  point,  then  scaling  by  1/2$^,  where  is  the  magnitude  of 
the  displacement  sensitivity  in  the  j  direction.  For  convenience,  the  strain  data  were 
decimated  to  the  £2  step  size  in  each  £3  plane  in  order  to  have  equal  resolution  in  both 
the  £1  and  £2  directions.  The  strain  images  were  then  median  filtered  with  a  5  x  5 
window,  resulting  in  a  slight  decrease  in  spatial  resolution.  These  strains  were  used  as 
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input  for  the  elasticity  reconstruction. 

The  3D  elasticity  reconstruction  was  performed  using  the  least-squared  error 
minimization  algorithm  discussed  in  Skovoroda  et  al  (1999),  with  a  second-order,  one¬ 
sided  finite  derivative  approximation  in  the  x 3  direction.  The  reconstruction  of  /i(r)  is 
a  boundary  value  problem,  therefore  a  unique  solution  is  obtained  only  with  boundary 
conditions.  So  a  square  35-mm  by  35-mm  region  of  interest,  which  contained  the 
inclusion  in  several  X3  planes,  was  identified  in  the  aq  and  x2  directions.  Along  the 
boundaries  of  these  regions,  and  in  the  first  two  X3  planes  furthest  from  the  center  of 
the  inclusion,  the  value  of  the  shear  modulus  was  set  to  one,  resulting  in  a  relative  shear 
modulus  reconstruction. 

5.  Results 

Representative  magnitude  and  corrected  phase  images  of  the  Semicosil  phantom  for  a 
1.5-mm  thick  plane  centered  about  X3  =  0.75  mm,  or  X3/R  ~  0.05,  are  shown  in  figure  3. 
Knowing  that  «  5.11  n  mm-1  in  the  aq  and  x2  directions,  the  number  of  2n  phase 
wraps  in  figure  3(c)  indicates  a  vertical  deformation  of  approximately  2.3  mm,  and  those 
in  figure  3(b)  a  horizontal  deformation  of  about  2.0  mm.  Reduced  phase  slope  in  the 
region  of  the  hard  inclusion  is  clearly  visible  in  these  figures  as  well.  Due  to  the  central 
location  of  this  plane,  there  is  little  feature  in  $3  (part  (d)). 

Figure  4  shows  representative  strain  maps  from  the  planes  centered  around  x3  = 
15.75  mm  and  x3  =  14.25  mm.  Due  to  the  loaded  state  of  the  phantom  during  imaging, 
the  sphere  became  prolate,  therefore  these  planes  correspond  to  the  X3/R  fa  1.05  and 
X3/R  fa  0.95  locations,  respectively.  One  normal  strain,  £22  (parts  (a)  and  (e)),  the 
in-plane  shear  strain,  £12  =  £21  (parts  (b)  and  (f)),  and  one  through-plane  shear  strain, 
£13  =  e31  (parts  (c)  and  (g)),  are  shown  for  each  plane.  These  components  are  all 
required  to  perform  the  elasticity  reconstruction  in  (8).  Note  that  the  presence  of 
through-plane  strains  in  (8)  necessitates  measurement  of  the  full  3D  displacement  field. 
In  addition,  although  elasticity  specific  details  are  seen  in  the  strain  maps,  features 
related  to  geometry  and  the  applied  deformation  are  also  clearly  present.  This  points 
to  the  need  for  a  proper  elasticity  reconstruction  in  order  to  obviate  these  confounding 
factors.  Also  shown  are  the  traces  of  the  strain  tensor,  £n  +  £22  +  £33,  for  both  planes 
(parts  (d)  and  (h)).  The  relative  lack  of  features  in  the  trace  of  the  strain  tensor  indicates 
that  the  phantom  is  nearly  incompressible  (like  soft  tissue). 

Magnitude  images  of  the  35-mm  by  35-mm  regions  of  interest  in  the  same  two 
planes,  along  with  two  different  shear  modulus  reconstructions  of  these  planes,  are 
presented  in  figure  5.  As  clearly  seen  in  the  magnitude  images,  the  hard  inclusion  is 
present  in  the  x3/R  fa  0.95  plane,  while  it  is  absent  in  the  x3/R  fa  1.05  plane.  Figures 
5(b)  and  (f)  show  2D  elasticity  reconstructions  of  these  two  planes,  while  figures  5(c) 
and  (g)  show  the  corresponding  3D  elasticity  reconstructions.  As  in  figure  1,  one  sees 
an  over-estimate  of  the  3D  spatial  extent  of  the  inclusion  in  the  2D  reconstructions. 
This  over-estimate  is  corrected  with  the  3D  reconstruction.  For  ease  of  comparison, 
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(a)S1(x1,x2)  (b)<t>1(x1,x2) 


Figure  3.  Representative  magnitude  and  phase  images  from  the  X3/R  &  0.05  plane 
of  the  3D  displacement  encoded  data  set  from  a  phantom  with  a  single,  hard,  spherical 
inclusion.  Si,  (a),  is  the  magnitude  of  the  a;  1 -displacement  encoded  data,  and  4>i, 
(b)-(d),  are  the  phase  images  of  the  2, -displacement  encoded  data. 


vertical  profiles  through  the  center  of  the  inclusion  from  the  magnitude,  2D,  and  3D 
reconstructions  are  presented  in  figures  5(d)  and  (h).  Note  that  the  magnitude  images 
only  convey  geometric  information. 

6.  Discussion 

The  stimulated  echo  sequence  presented  phase  encodes  internal  displacements  using 
gradient  pulses.  An  externally  applied  deformation,  synchronized  with  the  pulse 
sequence,  produces  an  internal  displacement  field.  This  deformation  is  actively  driven 
with  a  pneumatic  device,  and  the  mechanical  transition  from  pre-  to  post-deformation 
occurs  during  the  sequence  mixing  time,  Tm-,  while  the  relevant  magnetization  is 
longitudinal.  Because  longitudinal  magnetization  decays  only  as  7\,  the  mechanical 
transition  period  may  be  extended  to  allow  potentially  long-lived  or  ill-defined  motions 
within  the  object  to  dampen.  With  a  sufficiently  long  Tm-,  the  encoded  displacement 
will  be  approximately  static.  However,  signal  loss  due  to  T\  relaxation  sets  a  practical 
limit  on  the  length  of  Tm-  To  determine  an  appropriate  mixing  time,  a  series  of  2D 
displacement  encoded  images  was  taken,  varying  Tm  from  50-750  ms.  From  the  phase 
maps  of  these  data,  one  could  see  that  the  top  acrylic  plate  of  the  deformation  device 
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(e)  e22(xrx2) 


(f)  e12(xvx2) 


(9)  e13(x1,x2) 


(h)  e11+e22+e33 


Figure  4.  Representative  strain  images  from  the  xz/R  ps  1.05,  (a)-(d),  and 
X3/R  rs  0.95,  (e)-(h),  planes  from  the  3D  displacement  encoded  data  of  the  phantom 
with  a  single,  hard,  spherical  inclusion.  One  longitudinal  strain,  £22,  the  in-plane  shear 
strain,  £12  =  £21  >  one  through-plane  shear  strain,  £13  =  £31,  and  the  trace  of  the  strain 
tensor,  £n  +  £22  +£33,  are  presented  for  each  plane.  Linear  scales  for  each  image  are, 
from  black  to  white:  (a),  (e):  [-0.08,  0.03];  (b),  (f):  [-0.05,  0.05];  (c),  (d),  (g),  (h): 
[-0.03,  0.03]. 


completed  its  excursion  in  under  200  ms,  and  the  internal  motion  in  the  phantom  became 
negligible  by  250  ms.  To  ensure  that  (3)  was  reasonably  satisfied,  a  270  ms  mixing  time 
was  chosen  for  subsequent  data  collection.  Although  the  deformation  device  used  here 
provides  adequate  transition  speed,  an  even  faster  device  would  allow  a  shorter  Tm, 
yielding  more  signal.  While  absent  in  the  phantom  used  here,  water  diffusion  in  the 
presence  of  displacement  encoding  gradients  will  be  another  source  of  signal  loss  in  in 
vivo  experiments.  This  loss  can  be  mitigated  by  reducing  the  displacement  sensitivity, 
$d,  or  by  shortening  Tm- 

The  quality  of  the  shear  modulus  reconstruction  ultimately  depends  on  the  local 
phase,  as  in  (11),  induced  by  the  encoding  gradients  and  the  local  displacement.  More 
specifically,  the  quality  depends  on  the  spatial  derivatives  of  the  encoded  phase.  A 
study  of  the  effects  of  displacement  sensitivity,  applied  deformation,  relative  hardness, 
and  diffusion  loss  on  the  signal-to-noise  ratio  (SNR)  of  the  elasticity  reconstruction 
has  been  presented  in  (Steele  et  al  1999).  This  study  demonstrates  that  increased 
intra-voxel  phase  wrap  will  increase  the  reconstruction  SNR,  up  to  a  tt  intra-voxel 
phase  distribution.  Note  that  the  reconstruction  SNR  increases  despite  a  reduction 
in  the  nuclear  magnetic  resonance  (NMR)  signal  from  the  object.  Assuming  a  linear 
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(a)s1(x1>x2)  (b)  k2(x1  ,x2)  (c)  k3(x1  ,x2)  (d)  Central  profiles 


x,  10  20  30 

x2  [mm] 


(e)S1(x1>x2) 


(f)  ic2(xi  ,x2)  (9)  k3(xvx2)  (h)  Central  profiles 


10  20  30 

x2  [mm] 


Figure  5.  an-displacement  encoded  magnitude  images,  Si  (^1,^2),  and  corresponding 
2D,  ^2(^1,  £2),  and  3D,  kz(xi,  £2)  elasticity  reconstructions  from  the  X3/R  «  1.05, 
(a)-(d),  and  X3/R  &  0.95,  (e)-(h),  planes  of  a  phantom  with  a  single  hard,  spherical 
inclusion.  Also  presented  are  the  central  vertical  profiles  of  each  distribution,  where 

( . )  is  «2j  and  ( - )  is  K3.  All  are  presented  on  a  log  scale  where  black  corresponds 

to  a  relative  shear  modulus  of  0.5  and  white  to  4.5.  The  background  has  a  realtive  shear 
modulus  of  1,  and  the  inclusion,  4,  from  independent  measurements.  For  geometric 
reference  purposes,  parts  (d)  and  (h)  include  plots  of  S^1  (filtered  and  normalized)  as 

( - )• 


phase  distribution  of  0  radians  across  a  voxel,  the  signal  modulation  from  that  voxel 
will  be  |sinc(0/2)|  =  | sin(0/2)(0/2)_1|.  However,  the  phase  gradients  (that  is,  the 
displacement  derivatives)  will  be  maximized  without  aliasing  as  the  intra-voxel  phase 
wrap  approaches  7r,  and  this  is  the  signal  that  is  important  in  the  reconstruction.  A 
7r  phase  wrap  may  be  achieved  through  many  combinations  of  applied  deformation 
and  displacement  sensitivity.  However,  increasing  will  increase  signal  loss  due  to 
diffusion,  as  discussed  above.  Hence,  a  smaller  displacement  sensitivity  and  increased 
deformation  would  appear  to  be  optimal.  Again,  there  is  a  trade-off:  as  deformation 
increases,  the  model  of  linear  elasticity  discussed  in  section  2.1  will  become  less  and  less 
valid.  Elasticity  reconstructions  from  finite  displacement  fields  have  been  demonstrated 
in  Skovoroda  et  al  (1999),  but  these  are  obviously  more  computationally  intensive  than 
the  linear  reconstructions  used  here.  In  relation  to  the  data  presented  here,  the  number 
of  2n  phase  bands  across  the  phantom  in  figure  3  clearly  indicate  that  this  data  was 
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acquired  with  a  sub-optimal  displacement  sensitivity /applied  deformation  combination. 
Because  neither  the  encoding  nor  the  deformation  used  here  were  extreme,  the  elasticity 
reconstruction’s  SNR  should  be  relatively  easily  improved  merely  by  optimizing  the 
intra- voxel  phase  wrap. 

Relative  hardness,  object  geometry,  and  deformation  geometry  also  effect  the 
displacement  phase  gradients.  In  general,  the  phase  gradients  increase  near  soft/hard 
interfaces  and  are  higher  in  relatively  soft  regions  of  tissue.  Excessive  phase  wrap 
(i.e.  strain)  can  lead  to  regions  of  significant  signal  loss  in  a  manner  analogous  to  flow 
dephasing  in  conventional  MRI.  The  resulting  reconstructions  would  suffer  from  this 
signal  loss.  Hence,  the  applied  deformation  and  displacement  sensitivity  should  be 
optimized  for  the  regions  of  highest  strain  in  an  object.  Increased  intra-voxel  phase 
wrap  in  regions  of  lower  strain  may  be  obtained  by  integrating  the  signal  from  several 
voxels;  in  essence,  applying  an  adaptive  voxel  size  based  upon  local  phase  gradients. 
This  increase  in  signal  would  come  at  the  expense  of  spatial  resolution.  Additionally, 
signal  loss  in  the  displacement  encoded  magnitude  images  may  be  useful  for  identifying 
high  strain  regions  in  tissue. 

Another  factor  affecting  the  displacement  signal  is  the  reproducibility  of  the 
applied  deformation.  For  multi-step  acquisitions,  as  those  presented  here,  good 
deformation  reproducibility  is  essential.  Variations  in  the  applied  deformation  will 
lead  to  phase  instability,  motion-like  artifacts,  and  errors  that  will  propagate  through 
the  elasticity  reconstruction.  Adequate  reproducibility  has  been  achieved  with  the 
current  deformation  system.  However,  irreproducible  or  asynchronous  motions  within 
the  imaged  object  may  be  problematic.  These  would  include  physiologic  cardiac 
and  respiratory  motion  present  in  in  vivo  experiments.  In  some  ways,  the  problems 
associated  with  undesired  motion  would  be  similar  to  those  encountered  in  diffusion 
MRI.  Because  the  applied  deformation  is  external,  though,  the  displacement  encoding 
can  be  tailored  to  it,  reducing  the  effect  of  undesired  motion  on  the  displacement 
data.  Additional  complications  arise  because  phase  derivatives  of  the  displacement 
data,  approximated  by  finite  differences,  are  required  for  reconstruction.  In  addition  to 
choosing  an  appropriate  deformation/encoding  combination,  methods  should  be  devised 
to  reduce  the  effects  of  undesired  motion  on  the  displacement  derivatives. 

Clearly  several  advantages  justify  performing  a  3D  elasticity  reconstruction  rather 
than  a  2D  reconstruction.  As  illustrated  in  figures  1  and  5,  and  as  discussed  in  sections 
2.2  and  5,  a  3D  reconstruction  provides  a  more  accurate  representation  of  the  elasticity 
distribution  than  a  2D  reconstruction  in  the  simple  phantom  used  here.  Complicated 
in  vivo  geometries  will  only  increase  the  likelihood  that  neglecting  out-of-plane  strain 
components  will  result  in  an  inaccurate  elasticity  estimate.  This  increased  accuracy 
comes,  though,  at  the  expense  of  increased  computational  complexity  and  increased 
scan  time.  This  increase  in  scan  time  may  be  lessened  through  the  use  of  echo-planar 
imaging  (Mansfield  and  Pykett  1978)  or  fast  spin-echoes  (Hennig  et  al  1986)  for  spatial 
encoding  (Chenevert  et  al  1999).  It  should  also  be  noted  that  the  reconstruction  in  (8) 
does  not  rely  on  the  assumption  of  incompressibility,  although  making  this  assumption 
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provides  another  means  of  regularizing  the  reconstruction. 

Additionally,  a  reconstruction  of  static  displacement  data  offers  several  advantages 
over  a  reconstruction  of  dynamic  displacement  data.  A  static  reconstruction  allows 
one  to  ignore  visco-elastic  effects  as  well  as  the  longitudinal  or  shear  nature  of  the 
applied  deformation.  Static  methods  also  provide  high  SNR  displacement  and  strain 
estimates.  Dynamic  methods,  on  the  other  hand,  provide  a  potentially  very  simple 
reconstruction  (Muthupillai  et  al  1995).  However,  as  previously  noted  in  section  1.3.1, 
this  reconstruction  may  be  compromised  by  interference  from  elastic  inhomogeneities, 
attenuation  of  shear  waves,  mixing  of  longitudinal  and  shear  waves,  and  resolution  limits 
imposed  by  the  shear-wave  wavelength.  Reconstruction  models  that  include  visco-elastic 
effects  allow  a  more  accurate  interpretation  of  dynamic  data  (Sinkus  et  al  1999),  (Van 
Houten  et  al  1999),  but  these  are  necessarily  more  complicated  than  static  models 
(Skovoroda  et  al  1995a,  1999). 

Choosing  a  contour  of  constant  shear  modulus  for  appropriate  boundary  conditions 
for  (8),  though,  can,  in  practice,  be  a  challenge.  In  the  applications  discussed  here, 
a  priori  knowledge  of  phantom  geometry  was  employed  in  the  reconstructions.  This 
may  be  possible  in  vivo  as  well,  albeit  more  complicated.  For  instance,  in  breast 
elasticity  imaging,  such  a  contour  may  be  defined  in  the  subcutaneous  fat  surrounding 
the  parenchyma  using  the  boundary  detection  procedure  described  in  Skovoroda  et  al 
(1995a).  The  elasticity  reconstruction  would  then  be  relative  to  the  shear  modulus  of 
the  fat  boundary,  assuming  that  it  is  constant.  Alternatively,  a  high  signal  cuff  of  known 
elastic  modulus  could  be  used  to  surround  the  breast.  This  would  provide  an  absolute 
image  of  shear  modulus  variations  if  the  boundary  contour  were  chosen  inside  the  cuff. 

The  3D  shear  elasticity  reconstructions  presented  above  contain  artifacts  both  inside 
and  outside  of  the  inclusion  due  to  the  finite  SNR  in  the  measured  displacement  strain 
components,  and  due  to  the  step  size  used  in  the  finite  approximation  to  the  derivatives 
in  (8).  In  contrast  to  the  2D  elasticity  reconstruction,  where  the  elasticity  distribution 
is  reconstructed  independently  in  each  plane,  the  3D  reconstruction  uses  the  elasticity 
distribution  in  neighboring  planes.  Therefore,  in  addition  to  in-plane  error  propagation 
problems  discussed  elsewhere  (Skovoroda  et  al  1995a),  error  propagation  in  the  through- 
plane  direction  may  occur  due  to  inaccurate  elasticity  reconstructions  in  the  preceding 
planes.  This  is  particularly  true  if  the  3D  elasticity  reconstruction  is  performed,  as  in 
this  paper,  by  solving  an  initial  value  problem  in  the  through-plane  direction.  Even 
though  the  more  stable  integral  based  approach  (Skovoroda  et  al  1999)  was  employed 
to  solve  (8)  for  each  plane,  the  results  of  the  3D  elasticity  reconstructions  in  subsequent 
planes  exhibit  significant  error  propagation  in  the  x3  direction. 

Given  a  particular  spatial  discretization  of  the  displacement  data,  the  error 
propagation  can  be  reduced  by  several  approaches.  These  include  more  appropriate 
data  filtering  and  choice  of  boundary  conditions,  but  these  considerations  are  beyond 
the  scope  of  this  paper. 

These  include  more  appropriate  data  filtering  and  development  of  less  sensitive  to 
noise  elasticity  reconstruction,  but  these  considerations  are  beyond  the  scope  of  this 
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paper. 

7.  Conclusions 

The  ultimate  goal  of  quantitative  elasticity  imaging  is  to  provide  physicians  with  a 
method  of  remotely  palpating  soft  tissue  to  detect  disease.  The  three  dimensional 
elasticity  imaging  technique  demonstrated  here  is  a  step  toward  extending  the  range 
and  sensitivity  of  palpation,  a  powerful  diagnostic  tool.  One  possible  application  of  this 
technique  would  be  measuring  the  elasticity  of  breast  tissue  normally  inaccessible  to 
manual  palpation.  A  large  elastic  modulus  difference  between  normal  and  pathological 
breast  tissue  has  been  measured  in  situ.  A  previously  mentioned  study  indicates  that 
soft  tissues  in  different  physiologic  states  display  shear  modulus  variations  of  1-2  orders 
of  magnitude  (Sarvazyan  et  al  1995).  If  these  elastic  changes  predate  calcification 
formation,  elasticity  imaging  may  increase  sensitivity  to  and  characterization  of 
malignant  breast  masses,  complementing  existing  diagnostic  tools.  The  relatively  high 
cost  of  MRI  may  hinder  using  this  approach  as  a  general  screening  technique.  However, 
additional  work  to  define  the  role  of  this  modality  as  a  primary  or  complementary 
diagnostic  tool  in  diseases  of  soft  tissues  seems  worthwhile  indeed. 
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INTRODUCTION 

Elasticity  imaging  may  offer  clinicians  a  non-invasive 
method  of  remotely  palpating  patients  to  detect  diseased  tissues. 
Tissue  elasticity  information  may  be  obtained  via  NMR  detection 
of  shear  waves1  or  of  static  tissue  displacement2  Here  we 
present  simulation  results  elucidating  the  effects  of  displacement 
sensitivity  and  applied  deformation  on  the  signal-to-noise  ratio 
(SNR)  of  elasticity  images  obtained  using  static  displacement, 
stimulated  echo  NMRI  (SDSEI). 

PRINCIPLE 

In  SDSEI,  the  object  is  deformed  during  the  stimulated  echo 
mixing  time,  TM.  Local  displacements  are  encoded  by  means  of 
a  pulsed  gradient  applied  before  and  after  deformation.  The 
sensitivity  is  given  by 

®d=YGdT 

where  the  displacement  sensitivity,  Od,  is  in  radians/distance,  Gd 
is  the  displacement  encoding  gradient  strength,  and  ris  the 
displacement  encoding  gradient  duration.  The  gradient  of  the 
displacement  field  yields  the  strain  field,  and  the  strain  field  is 
used  as  an  input  for  a  boundary  value  problem  to  extract  the 
tissue  elastic  moduli.3 
SIMULATION 

A  simple  one-dimensional,  100-mm  long  object  is  used  in 
all  simulations.  This  object  contains  an  inclusion  whose  shear 
modulus  is  five  times  larger  than  the  surrounding  material’s. 

The  pulse  sequence  assumed  is  that  in  the  paper  by  Chenevert  et 
al. ,2  with  TM  =  200  ms,  x  =  4.5  ms,  and  225  ms  between 
displacement  encoding  gradients  .  The  NMR  resolution  is  128 
pixels.  The  base  SNR  of  the  material  is  set  to  130,  matching 
the  SNR  of  a  tissue-mimicking  material  measured  using  the 
SDSEI  pulse  sequence. 

The  variables  used  are  Gd  and  the  normalized  applied  surface 
deformation,  Gq.  Since  x  is  assumed  constant,  <£d  varies  linearly 
with  Gd.  Gd  ranges  from  0  to  5  G/cm  in  0.2  G/cm  increments, 
and  Go  from  1  %  to  20%  in  1%  steps.  For  each  gradient- 
deformation  combination,  both  the  NMR  image  and  the  elasticity 
image  are  reconstructed.  The  reconstructed  elastic  moduli  are 
normalized  to  the  background  reconstruction.  Then  the  NMR 
image  SNR,  SNRo,  and  the  elasticity  image  SNR,  SNRE,  are 
calculated.  Four  different  noise  instances  are  averaged  for  every 
point. 

RESULTS 

Figures  la  and  lb  show  the  SNRo  and  the  SNRE, 
respectively,  of  the  background.  Here  we  see  that  the  SNRo 
decreases  with  increasing  Gd  and  Go  due  to  intrapixel  dephasing. 
Note  that  signal  loss  due  to  diffusion  increases  as  Gd  increases. 
The  SNRE  begins  at  0  for  Gd  =  0,  and  then  increases  with  both 
deformation  and  sensitivity  until  the  intrapixel  dephasing  is  n 
radians.  Beyond  this  point  it  drops  off  quickly  to  zero.  Again, 
diffusion  signal  loss  is  present 

Figures  2a  and  2b  show,  respectively,  the  SNRo  and  the 
SNRe  of  the  inclusion.  Here  the  major  source  of  signal  loss  in 
SNRo  is  diffusion.  Because  the  inclusion  is  harder  than  the 
background,  it  tends  to  displace  rather  than  deform  under  the 
applied  surface  deformation;  hence  only  the  gradient  contributes 
significantly  to  the  intrapixel  dephasing.  The  SNRE  simply 
seems  to  exhibit  diffusion  loss  and  does  not  peak  as  the  SNRE  of 
the  background  does.  Also  note  that  the  SNRE  of  the  inclusion 
is"  lower  than  that  of  the  background. 
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This  can  be  explained  by  noting  that  the  elasticity 
reconstruction  depends  upon  the  strain  field,  and  the  strain  is  the 
gradient  of  the  displacement  field.  As  noted  above,  the  hard 
inclusion  tends  to  displace,  not  deform,  under  a  certain  surface 
deformation.  Thus,  not  only  does  the  intrapixel  dephasing 
decrease,  but  the  m/e/pixel  dephasing  (i.e.,  strain)  decreases. 
Hence  the  combinations  of  Od  and  e0  presented  here  do  not 
approach  the  optimal  strain  fields  for  the  inclusion 
reconstruction. 


Fig.  la:  Background  SNRo 
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Fig.  2a:  Inclusion  SNR0 
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Fig.  2b:  Inclusion  SNRE 
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Encoding  gradient  (G/cm) 

CONCLUSIONS 

We  have  illustrated  the  effect  of  displacement  sensitivity  and 
applied  deformation  on  static  displacement  NMR  elasticity 
images  and  noted  the  differing  optimal  imaging  parameters  for 
tissues  with  different  shear  moduli.  These  results  suggest  a  type 
of  adaptive  elasticity  reconstruction  wherein  the  voxel  size  for  a 
given  region  of  tissue  varies  inversely  with  the  strain  field  in 
order  to  optimize  the  elasticity  reconstruction’s  SNR. 
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